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ABSTPIA.CT 

To aid in the modeling of a steady state spin, the 
equations of motion of an airplane are formulated in a 
cylindrical coordinate reference frame. The derivation 
of the equations is presented and the resulting equa- 
tions are simplified for the equilibrium spin condition. 
These simplified equations are used in an unconstrained 
computer parameter optimization technique that algebrai- 
cally solves the differential equations for the equilib- 
rium state. The results of the computer work are 
presented and compared with previous prediction schemes. 
The potential of the method is demonstrated by applica- 
tion to a study of the effects of density variation. 
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LIST OF SYMBOLS 



The definitions 


of the symbols used in this paper are as 



follows : 



b 


Wing span, ft 


c 


Mean aerodynamic chord, ft 


F 


Lagrange generalized force, lbs or 
ft-lbs 




Aerodynamic forces along the body axes, 
lbs 


g 


2 

Acceleration of gravity, ft/sec 


^X'^Y'^Z'^XZ 


Body axes moments and products of 
inertia, slugs— ff^ 


M 

X 


Rolling moment acting about the X 
body axis, ft— lbs 


My 


Pitching moment acting about the Y 
body axis, ft-lbs 


^z 


Yawing moment acting about the Z 
body axis, ft-lbs 


m 


Airplane mass, slugs 


P/q^r 


Angular rates of rotation about the 
X, Y and Z body axes, rad/sec 


^cx, 


2 

Dynamic pressure, 1/2 pV 


^r 


Lagrange generalized coordinate 


R 


Radius of helical path of aircraft, ft 


S 


2 

Wing area, ft 


T 


Lagrange kinetic energy, ft-lbs 


V 


Inertial velocity of the C.G., ft/sec 


V 


Lagrange potential energy, ft— lbs 


X,Y,Z 


Aircraft body axes 
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Cartesian coordinate reference axes 
that are attached to the cylindrical 
radius vector. 

Aircraft altitude, ft 

Angle of attack, degrees 

Direction cosine, dimensionless 

Angle of sideslip, degrees 

Cylindrical orientation angle, radians 

Aileron deflection (positive when right 
aileron trailing edge is down) , degrees 

Elevator deflection (positive when 
elevator trailing edge is down) , degrees 

Rudder deflection (positive when rudder 
trailing edge is left when viewed from 
above) , degrees 
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Atmospheric density, slugs/ft 
Euler orientation angles 
Derivative with respect to time 
and derivatives: 
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INTRODUCTION 



I . 



A. BACKGROUND 

As a result of post— stall/spin accidents a significant 
number of lost fighter aircraft and pilot fatalities have 
been recorded over the last decade. Investigations have 
shown that these aircraft exhibit poor spin characteristics 
and that recovery from a fully— developed spin is usually 
difficult or impossible. For these reasons the spin is no 
longer considered a tactical maneuver and is regarded by 
pilots as an undesirable and potentially dangerous flight 
condition . 

Although the flight conditions leading to a spin can 
generally be avoided, the combat pilot often finds the need 
to utilize the edge of the flight envelope in a combat 
situation. More often than not, it is only in, this boundary 
area that the combat superiority of one aircraft over an- 
other is realized. If the aircraft demonstrates unsatis- 
factory stall and spin characteristics, the pilot is 
disinclined or restricted from utilizing this outermost 
region of flight to his best advantage. As a result, the 
overall weapon system effectiveness and margin of superi- 
ority of the aircraft is lost. 

In light of the poor stall and spin characteristics 
of current fighter configurations a more detailed consid- 
eration must be given to these characteristics during the 



10 



early design stages of future fighter aircraft. In recog- 
nition of this fact this research was undertaken to provide 
a more accurate analytic tool with which the designer can 
investigate the spin characteris tics of future aircraft and 
to aid in understanding how various factors in design may 
affect these spin motions. It was not proposed that this 
effort solve the spin problem itself, but it was hoped that 
it would prove to be a useful approach and that it would 
provide information complementary to that obtained experi- 
mentally . 

B. SCOPE 

Traditionally, analytic spin prediction has involved 
either an approximated solution (by requiring that only 
part of the equilibrium conditions be satisfied) or by com- 
puter generation of time— histories of characteristic para- 
meters (by integrating the equations of motion forward in 
time until the average values of the variables are approxi- 
mately constant) . More recently the idea of mathematically 
casting the equations of motion into a form specifically 
intended for spin prediction has shown some promise. Of 
particular note is the work carried on by Buehler (Ref. 1) 
and Champoux (Ref. 2). These earlier works provided much 
of the initial groundwork for this research. 

In developing the analytic approach to the spin problem 
a set of objectives was established to subdivide the re- 
search problem and also to serve as guidelines to ensure 
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that the work proceeded in the proper direction as effec- 
tively as possible. These objectives were: 

1. To suitably model the aircraft in a spin. 

2. To develop the spin equations of motion. 

3. To establish a suitable technique of solving the 
equations of motion for the equilibrium spin con- 
dition . 

4. To develop a computer program that will predict 
aircraft spin characteristics. In addition, the 
program should be flexible enough to aid in the 
understanding of how various parameters effect the 
developed spin motions. 

5. To demonstrate the usefulness of the method and to 
verify the results with other prediction schemes. 
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II. DEVELOPMENT OF THE AIRCRAFT EQUATIONS OF MOTION 



A. COORDINATE SYSTEM 

To an observer on earth, the aircraft's center of grav- 
ity in a steady spin appears to trace a descending constant 
radius helical path about an imaginary vertical central 
axis. By assuming that the rotation of the earth is negli- 
gible, an inertial reference frame can be established that 
is fixed to the earth's surface at the point of intersection 
of this imaginary spin axis. 

The most convenient method of describing this motion is 
through the use of cylindrical coordinates. In this formu- 
lation, the central axis of the cylinder is superimposed 
on the imaginary aircraft spin axis as shown in Figure 1. 

The advantage of this choice for an inertial coordinate sys- 
tem is that the center of gravity of the aircraft is con- 
veniently located in terms of an altitude coordinate (Z^) , 
a radial coordinate (R) and an angular position coordinate 
(Y). 

The body axes, as shown in Figure 2, were chosen for 
the formulation of the equations of motion. This allows 
the direct use of standard wind tunnel data which was avail- 
able- for this research effort. Symmetry of the aircraft 
about the X and Z axes was assumed so that all cross-product 

of inertia terms, with the exception of I , were identi— 

XZ 

cally zero. 
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Figure 1. Inertial Cylindrical Coordinate 
System. 




Figure 2. Aircraft Axes. 
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The orientation of the aircraft is described by intro- 
ducing another coordinate system and a set of three rota- 
tion angles. A fixed cartesian coordinate reference frame 
(Xi/YijZf) is positioned at the center of gravity of the 
aircraft and aligned with respect to the inertial cylindri- 
cal frame so that the axis is parallel with the vertical 
central axis and the axis coincides with the — R position 
vector. Although this cartesian coordinate reference frame 
remains fixed in orientation with respect to the inertial 
cylindrical frame, it is still free to move in space as Z^, 

R and y take on different values. The orientation of the 
body axes relative to the cylindrical frame of reference is 
given by the aircraft Euler angles (^,6,4>), which define an 
ordered series of three consecutive rotations as shov/n in 
Figure 3. The three position coordinates (Z^,R,y)/ together 
with the three orientation angles allow a complete 

description of the motion of the aircraft in space and will 
serve as the generalized coordinates in the Lagrangian 
development of the aircraft equations of motion. 

B. AXES TRANSFORMATION AND AIRCRAFT ANGULAR RATES 

Before proceeding with the development of the equations 
of motion, certain algebraic relations are needed that will 
transform the actual forces, moments, velocities and angular 
rates of the aircraft, about its body axes, to the cylindri- 
cal reference frame. These transformations are nothing more 
than vector projections from one system of axes onto another. 



15 




Figure 3. Spin Model. 

Consider that the body axes are initially coincident 
with the inertial frame and then are rotated 

about by an angle if) to form an intermediate reference 

p 

frame , as shown in Figure 4 

Since both of these orthogonal co- 
ordinate systems form a basis for 
three— space, any point (P) can be 
described by a vector (R) from the 
origin to that point in terms of 
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either system. 

Mathematically, this can be expressed: 



r ^2^2 ^2^2 ^2^2 



( 1 ) 
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( 2 ) 


k2‘ 


r = Z 2 


= x^k^- 


/N 

^1 




Y k • 
12 


^1 




"A- 


/N 

^1 





which is written more compactly in matrix form as: 



/S /\ /S yN /\ 




The three— by— three matrix of unit vectors is called a trans- 
formation matrix and denoted, in this case, by t'^2/1^ 
read "the transformation from frame one to two." Carrying 
out the indicated dot products for the rotation, the trans- 
formation matrix has the value: 




cos ip —sinip 0 

sin \p cos ip 0 

0 0 1 



(4) 



The second rotation, 0, about the Y 2 axis carries sub- 
scripted coordinate system two into system three, as 
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shown in Figure 3. The resulting transformation matrix 
is : 




cos 6 
0 

—sin 6 



0 

1 

0 



sin 6* 
0 

cos 9 



(5) 



The final rotation, <(> , about the axis carries sub- 
scripted coordinate system three into the unsubscripted 
body axes, as shown in Figure 3. This resulting trans- 
formation matrix is: 




"1 

0 

0 



0 0 - 

cos (|) sin ({) 

, — sin 4 > cos (j) 



( 6 ) 



These individual transformation matrices can be multiplied 
together to form a single three— by— three matrix which will 
be denoted by [T ] • Since each transformation matrix 

p/1 

premultiplies the vector arrived at in the previous step, 
the total transformation from the cartesian coordinate 
frame to the aircraft body axes is given by: 

[Vi] “ ["p/3] 



3/2 



2/1 



(7) 
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Carrying out the matrix multiplication, the total trans- 



formation has the following form in which denotes its 

elements : 



‘Vi> “ 



a 



COS0COSlf) 


tti 2 = cos0sinif> 


— sin0 


„ = sin<f)sin9 
costfj 


a„ 2 “ sin(f)sinO 
sinifj 


^ 22 ~ sin(f)cos0 


— cos(f)Sim|^ 


+COS(f)COSlfi 




- = cos(f)Sin9 

COSlf) 


sinif) 


cos(f)cos0 



+sin(|)sin4» 



( 8 ) 



— sin(i)COsii> 



In relating the angular velocity rates of the aircraft, 
p, q and r to the inertial frame, we make use of the trans- 
formations ' ^'^p/3^ ^"^3/2^’ This development 

entails projecting the orientation angular rates, , 9 and 
<j) onto the aircraft body axes. In addition, the inertial 
angular rate (y) must be included to account for the rota- 
tion of the cartesian coordinate reference frame , Z^) . 

Both Y and are considered to be applied about the 
axis of the cartesian coordinate reference system and, 
through the transformation, their projections are 

realized. The Euler pitch rate (0), about the axis, in- 
volves only two transformations ^'^ 2 / 2 ^ ^'^p/3^' pro- 

jection on the body axes. The Euler roll rate ((f)) about 
the body axis involves no transformations. Mathematically 
expressed : 



( 9 ) 



"1 
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1 ° ' 




1 ^ 1 




0 




► + 


1 ° 1 


r I 


i if)+Y ‘ 


( ^ 




1 ® I 
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On substitution of Equations (5) , (6) and (8) for ["^'2/2^ ' 

respectively, in Equation (9), the exact 
relation between the aircraft angular rates and the Euler 
angle rates is given by: 

p = — sin6 (<i^+Y) 

q = 6 cos4> + COS0 sin(|) (ij;+Y)' (10) 

r = cos0cos(|) (ijj+Y) “0 sincj) 



The time rates of change of the direction cosines (*^j^j) 
and the angular rates (p,q and r) are obtained by direct 
differentiation and are presented in Appendix A. 



C. LAGRANGE FORMULATION OF THE AIRCRAFT EQUATIONS 
OF MOTION 

Assuming that the aircraft is a single rigid body, the 
three position coordinates (Z^,R,y) and the three orienta- 
tion angles (\p,6,<})) can completely describe the six degrees 
of freedom of the vehicle. These six variables form the 
independent set of generalized coordinates necessary for 
Lagrange modeling of the motion. Let T denote the kinetic 
energy of the system; V, the system potential energy; q^, 

a generalized coordinate; and F^ a generalized force due 

-ir 

to non— conservative forces. Then the Lagrange equations 
for this system can be written: 



d .3T 3T 

' 3q 
3q ^r 



•s 



r = 1 , 2 , ... 6 



( 11 ) 



The kinetic energy of the aircraft, relative to the 
inertial reference frame, is given by the superposition of 



20 



the energy of the center of gravity, plus inertial rota- 
tional terms about the center of gravity, and is written: 

T = |m[(YR)^+ (R)^+ (Zq)^] + 

It is necessary to use the relations given in Equation (10) , 
as Equation (12) is not compatible with the generalized co- 
ordinates (Z^,R,y / ^ f 9 , 4>) . Upon substitution of these rela- 
tions for p, q and r the compatible kinetic energy equation 
is ; 



T = ^[(yR)^+ (R)^+ (Zq)^] + 

1 * . * * 2 

+ coscj) + sin(j)cos0 (i(j+y) ] 



sinB (ip+Y ) ] 



(13) 



1 • • 2 

+ f‘^°^4>COS0 (l/j + Y ) ] 



+ [({) — sin0 (i|;+Y ) ] [(4 '+y) cos(})COs 0 — 0 sin<})] 

A/j 



If a reference for the potential energy is established 
at sea level, V is given by: 

V = niqZ^ (14) 

The partial derivatives of Equation (13) with respect to 

9 T 

the time rates of the generalized coordinates (— ;- ) are: 

9<ir 

— = mZ (15) 

3Z ° 

o 

3T 

9R 
3T 

3Y 



= mR (16) 

= mRY+ IxP“i3+ Iy^“23^ ^Z^^33 

+ I (p cos(})Cos6 — r sinO) (17) 

X ct 



21 



and 

9T 

3ii^ 

3T 

30 

3T 

3^ 



^xP“l3 ^y‘J“23 ■" ^Z'"“33 ~ ^XZP <“> 

I q cos(J) + I r sin<|) + (pcos(j)cos0— rsinG ) (19) 

Y Z X Z 

^xP + Ixz'! -- <2“' 



The time rates of change of Equations (15) through (20) 
^ d ,3T , 

represent (— ) and are: 

9qr 



d ,9T , " 

3Z^ ° 



( 21 ) 



d /9T. _ 

dt • 

3R 



(22) 



d 3T . . • . • 

^(— ) = mR(RY+ 2 RY)+I„(pa 3+pa )+I (qa23+qa23) + 

3y 

1 2 (ra33+ra33) +Ij^2 (~P (6cos(j)sin0+^sin(})cos0 ) + 



pcos(|)cos0— 0rcos0 — r sin0) 



(23) 



d 9 T 

dt<— > “ ^x'P“13'"P“i3>+^y‘'3“23+‘J“23>+^z'P“33+™33> + 



3i|; 



fv -7 (“P ( 0cos(j)sin0+(j)sin(j)cos0 ) +pcosc})cos0 — 
X z 



0 r COS0 — r sin0) 



(24) 



dt ' • ' 
30 



— (— ) 
dt ^ ' 

34) 



I (qcos4)— q^sin4> ) +I„ (rsin4)+r^cos4) ) 
^ z 

(P^cos4i+psin4)) 






(25) 

(26) 
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The partial derivatives of Equation (13) with respect 



9T 

to the generalized coordinates (-^ ) are: 

"r 



3T 

9Z 

o 



0 



3T 

9R 



• 2^ 
my R 



(27) 

(28) 



9T 

3Y 



= 0 



(29) 



9T 

dip 



= 0 



9T 

= —(ip+y) (I pcos6+I„qsin0sin<|)+I„rsin9cos<})) 
ou A Y Z 

• • • • 

+ I (— pcos(J)sin9 (i|j+y) — rcos9 (i|)+y ) ) 

9T 

= qrd^-I^) + “ sin4)COs9 (i|)+y) ) 



(30) 



(31) 

(32) 



The partial derivatives of Equation (14) with respect 

3 V 

to the generalized coordinates (-^ ) are: 

'^r 



9V 

o 

9V _ 9V _ 9V _ 9V _ 9V 

9R 9y 9i|; 99 94> 



(33) 

(34) 



Equations (21) through (34) complete the formulation 

of the left hand sides of the six Lagrange equations, as 

generally expressed by Equation (11) . In completing these 

equations it is only necessary to formulate the generalized 

forces F„ , F_, F , F,, F„ and F , . The principle of vir— 

Z^ R y i|j 9 (|) 

tual work forms the basis for these calculations and can be 
described mathematically by: 
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(35) 






My( 



9 e. 



9e, 



-) + M, ( 



Z'9q ^ 



where X,Y,Z are positive incremental displacements along 
the axes, e^,e^,e^ are positive incremental angu- 
lar displacements about the axes. 

The incremental translation or rotation along or about 
the 2^ axis, caused by varying each generalized 

coordinate individually, can easily be visualized with the 
aid of Figure 3. For example, varying the generalized co- 
ordinate Z^ by a positive increment will produce no trans- 
lation or rotation about either or and no rotation 

about Z^. Using Equation (35) as a guide (realizing that 
each X,Y and Z appearing should be subscripted with a one, 
because we are dealing with the cartesian coordinate refer— 



9X, 



9Y, 



"'■^1 ^^X ^^Y 

ence frame) , it follows that f and 

dZdZoZoZ 

de^ O O O O 

are identically zero. From Figure 5 it is clear that 
o 




Figure 5. Positive Incremental Displacement of the 

Generalized Coordinate Z . 

o 
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as is incremented a positive amount, follows direct- 
ly in its negative defined direction. This implies that 
= —1 and the generalized force F„ = — F„ follows from 



3Z. 



Equation (35) . 

The remaining generalized forces with respect to the 
cartesian coordinate reference frame (X^,Y^,Z^) can be 
derived in a similar manner. The complete set is: 

= 1 



F , = 

\p Z- 



Fq = — sin + M cos 
F, = M cosOcosi|; + M cosOsimp — M sin9 

(p Xj^ Y^ 



F = — F 
Z Z^ 

o 1 



F = — F 
R Y. 



F = F R + M 
Y Z^ 



(36) 

(37) 

(38) 

(39) 

(40) 

(41) 



Making use of the axes transformation , the general- 

ized forces expressed in terms of the body axes are: 



^R ^“l2^X °‘32^Y ^32^Z^ 



(42) 



^^°‘ll^x'*' “21^Y “31^Z^ ^“l3^X °‘23 ^y"'’ ^33^^Z^ 



F., = + ci23H^ + 



13 X 



'33 Z 



Fq = My cos (p — sin ip 



(44) 

(45) 



<p X 



^Z^ ^“l3^X ^23^Y “33^Z^ 



(46) 

(47) 
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Adding Equations (42) through (47) to the respective 
right-hand sides of the earlier mentioned grouping (of the 
left-hand terms) , in the general expression given by Equa- 
tion (11) , completes the development. This grouping is 
shown in Appendix A, where the six equations of motion are 
presented in their entirety-. 

Engine thrust terms are not included in the six equa- 
tions of motion. Previous studies indicate .that use of 
thrust in the spin is normally avoided since flameouts are 
likely and serious damage can result (Ref. 4). 

D. STEADY SPIN EQUATIONS 

The full equations of motion (as given in Appendix A) 
can be greatly simplified for a steady state spin condition. 
In such a state, the aircraft is imagined to be descending 
at a constant velocity, at a constant radius (R) and rota— 
tion rate (y) about the central spin axis and with a fixed 
orientation with respect to the cartesian coordinate refer- 
ence frame (Xj^ , , Z ^^) . This implies that all second— deriv- 

ative terms and any first— derivative terms of the Euler 
angles (\pf6,(p) can be set equal to zero. On examination 
of the auxiliary terms of Appendix A, it is easily seen that 
all time rates of the direction cosines angular 

accelerations (p, q, r) are equal to zero. After making 
these simplifications, the six equations of motion that de- 
scribe the steady spin condition are: 

RES(l) = 0 = + «23 Fy + 0133^2 ~ n^g (48) 
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RES (2) = 0 
RES (3) = 0 

RES (4) = 0 
RES (5) = 0 

RES ( 6 ) = 0 



The aerodynamic coefficients are written in a form con- 
sistent with the data obtained from Ref. 5. The aerodynam- 
ic forces and moments are v/ritten: 



^X = 


+ Sq • 'I' 

12, ® , 


(54) 
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JPV S(C^+C^ -6^+c^ -6^+c .p + -r) 

0 0 p r 

12 ^ ^ — 


(55) 
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JPV S(C + C -6^ + c -q) 

«e 

12, „ — — . 


(56) 
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ipV Sb(Cj^+ -6^+ Cj^-p + -r) 

• <5^ P r 

1 2 — — 


(57) 


My = 


4pV Sc C + C *6 + C ‘q) 

2^ m m. e m ^ 

12. . “ 


(58) 


^z = 


^pV Sb C +C *6+0 *6 +C ‘P+C T 

2"^ n n. a n. r n ^ n 

6 6 p r 

a r ^ 


(59) 



If the relations given by Equations (54) through (59) 
were substituted into the steady spin equations, as given 
by Equations (48) through (53), the six equations would 
then describe the resulting aircraft state in terms of con- 
ventional aerodynamic forces and moments. 



12 X 



“22^Y 


+ 


“32^Z - R 


(49) 


^X'^°‘21^Y 


+ 


“31^Z^ 




X ■*'^23^V 


+ 


“S3^Z^ 


(50) 



*^13^X *^23^V ^ 
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(51) 



— Y(IvPCOS0 + I qsinc()sin0 + 

X Y 

I (rcos<J)sin0) + M cos<}) — M sin(|) 

La ± /j 



(52) 



— yIv -7 (pcos(})sin0 + rcos0) 



M^-qr ( I 



(53) 
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The residual notation (RES(i)) introduced in the steady 
spin equations, as given by Equations (48) through (53), is 
used as a measure of how closely each associated equation 
is being satisfied. Ideally, each residual should be iden- 
tically equal to zero but this is not possible due to trun- 
cation and round-off error within the computer at various 
stages of the program. 

A criteria function is defined that will indicate how 
close the complete set of six equations is to a steady spin 
condition. The criteria function chosen for this investi- 
gation was the sum of the squares of the individual resid- 
uals, which is written: 

F = RES(1)^+ RES(2)^+ RES(3)^+ RES(4)^+ RES (5)^ + RES (6)^ 

(60) 

This function will be the object of the minimization effort 
in the optimization procedure described in the next section. 

The airplane has six degrees of freedom with second 
order differential equations; therefore, 12 variables are 

required to specify its state. For convenience the 12 

• • • • • • 

variables chosen are: ip , Jp , d , d , (p , <p , , Z^, R , R, y and y. Mak- 

ing use of the fact that each equality constraint imposed 
upon the system will remove one variable, the total number 
of required variables can therefore be reduced. Some of 
these have been discussed in describing the steady spin; 
specifically, that R=i|; = 0 = ^ = O. They were also used 
to reduce the general equations given in Appendix A to the 
steady spin equations given by Equations (48) through (53). 
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By specifying that the search be made at y = 0 and at a 
specified altitude (2^) / the total number of equality con- 
straints is six. This leaves six variables to search over 
if each control surface is held in a fixed position. The 
chosen set of variables to be searched over in this inves— 
tigation are: y, R, Z^, tp 9 3.nd (p . Specifying these six 

initial conditions, the altitude parameter (p, as 
doesn't appear explicitly in the steady spin equations), 
the physical parameters of the aircraft (mass, chord, span, 
etc.), and having a set of tabulated aerodynamic coeffi- 
cients, the criteria function can be evaluated. The only 
thing now required is a systematic method of varying the 
six variables (y, R, Z^, ip , 0, cj)) to obtain the minimum 
value of the criteria function close to zero. 

The solution of the steady spin condition involves the 
simultaneous solution of six highly coupled nonlinear dif- 
ferential equations involving aerodynamic data. There 
does not exist a closed form analytic solution to this set 
of equations and therefore a numerical method must be 
employed. 
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Ill . 



REVISED FORMULATION OF THE STEADY SPIN EQUATIONS 



The steady state equations were FORTRAN coded for use 
in program SPIN and initial efforts at predicting the equi- 
librium spin condition were unsuccessful. During these 
efforts it was realized that the rate of convergence was 
extremely slow and that the solutions were not close to the 
states already known to exist, as shown by Adams (Ref. 7) . 

After some analysis it was theorized that the deletion 
of the acceleration terms from the derived equations of 
motion was the cause. The equations should have been un- 
coupled and solved explicitly for these acceleration terms. 
Then the residual terms, representing the acceleration 
values, would be a true measure of just how close the air- 
craft is to the equilibrium spin condition. 

Attempts at uncoupling the full equations of motion, as 
given in Appendix A, showed that it would be, at most, a 
very time consuming and nearly impossible task as the ijj , 9 
and (p equations are so highly coupled. Instead, a search 
for an alternate set of moment equations was made with the 
decision to use the equations given by Etkin (Ref. 3) . The 
three equations express the moments about the body axes and 
are : 



M 



X 



= i^p 



- 






V - 



^XzP‘3 



(61) 
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( 62 ) 



”z =-^XzP + ^Z^ ■" P'3'^y“ ^Z> + ^xz'^r '“> 

The first three equations of the original formulation 
express Newton's second law of motion along the radial, 
tangential and vertical directions through the aircraft's 
center of gravity. These equations can be expressed, 
without introducing the transformations (a^j), about the 
cartesian coordinate reference frame as; 

F = my (64) 

^1 

F + my^R = mR (65) 

1 

F = mZ (66) 

^1 

These three force equations when combined with the three 
moment equations, as given by Equations (61) through (63), 
describe the motion of and rotation about the center of 
gravity of the aircraft and therefore completely describe 
the motion of the six— degree— of— freedom vehicle. 

By explicitly solving Equation (64) for the accelera- 
tion term y, Equation (65) for R, Equation (66) for Z and 
Equation (61) through (63) for p, q and r the equations 
are in the form desired for the optimization search. 
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IV. OPTIMIZATION PROCEDURE 



A. CONCEPT OF OPTIMIZATION 

Using state variable notation, the equilibrium solution 
and optimization scheme for’ this particular formulation can 
be described mathematically as follows. If we let 

S \ 

R 
Z 
P 

q 

r 

and 



X 



= < 






the acceleration terms 



therefore 
finds } Y 



= < 






X y = { F(Y) 

such that I 



, the independent 'variables , 



. Then the optimization scheme 
X I =0, for equilibrium. The op- 
timization scheme employed in this investigation uses the 
criteria function which is the sum of the squares of the 
acceleration terms (modified equation residuals), given by; 

J(Y) = X) x.^ = Y^+ R^+ Z^+ p^+ q^+ r^ 
i=l ^ 
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and finds | X. | such that J (Y) is an absolute minimum 
which yields | ^ | ~ 

B. OPTIMIZATION ALGORITHM 

The optimization method used in this investigation 
finds the minimum value of the criteria function, which is 
a real— valued function of several variables. The criteria 
function is used as an indicator of how close the aircraft 
is to being in an equilibrium spin condition. The minimi- 
zation process takes into account any constraints imposed 
upon the system, which, in this case, were imposed limits 
of angle of attack and sideslip from the available aero- 
dynamic coefficients. 

The optimization procedure starts with a given set of 
independent control parameters and varies the current values 
of these parameters, within the imposed boundaries, to im- 
prove the value of the criteria function. The process of 
parameter adjustment continues until the minimum of the 
criteria function is reached. Depending on the hyper— sur- 
face defined by the independent control parameters, the 
minimum found may be just a local minimum which would not 
satisfy the desired steady— spin criteria functional value 
(equal to zero) . The absolute minimum of the hyper— surf ace 
will yield a criteria function value equal to zero and is 
the object of the search. The initial conditions given 
for the independent control parameters will greatly influ- 
ence the success of finding the absolute minimum. There- 
fore, there is a certain trial and error process in the 
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search for the absolute minimum without any guarantee that 
one exists. This then is the main disadvantage of the op- 
timization procedure because in spite of the energy and 
perseverance of the researcher he cannot, at the end, say 
that a steady spin condition does not exist for that air- 
craft if in fact he hasn't been able to locate an absolute 
minimum. .This is a very common disadvantage shared by all 
optimization procedures available at this time. 

The optimization algorithm used in the spin program is 
called EXTREM. A complete listing of this subroutine is 
given in Appendix B along with a functional flow chart. The 
user specifies a set of independent variables in name and 
initial value and a corresponding number of scaled incre- 
ments defining the initial step sizes along each of the in- 
dependent variables. The algorithm first checks that the 
given initial conditions do not violate any constraints 
imposed on the system and then proceeds to determine the 
main line of search in the variable space. The initial con- 
ditions form a point in this variable space and the speci- 
fied initial step sizes give a second. A vector joining 
these two points determines the initial main line of search 
for the first stage, as shown in Figure 6 for a two— variable 
problem. Along this main line the approximate minimum point 
is found by extrapolating a parabola through three points 
determined by DX, X^^ and X^^+ DX and then interpolating 

for the minimum ^ Gram— Schmidt orthogonalization- 

process determines the secondary direction of the firsc 
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stage which is orthogonal to the main direction at ^ ' 

The approximate minimum point found along this line. 



XCx) 



iEC.C»lOW2.S t>lP.£CAI6M 




Figure 6. Main Line of Search in a Two Variable Problem 



For a third variable (not shown in Figure 6) a line of 
search through is determined by the Gram— Schmidt pro- 

cess that is orthogonal to both the main and secondary 
directions. Each new direction is orthogonal to all the 
preceding directions of search and the procedure for deter- 
mining these lines of search and approximate minimum points 

t h 

is continued until the minimum along the line (K being 

the number of independent variables) is found. At this 

point the first stage is complete and the main direction 

for the second stage is determined by a line joining the 

point and v;here is the approximate minimum 

t h. 

along the K line of the first stage. The second stage 
proceeds just as described for the first, and the overall 
stage procedure continues until a local or absolute rela- 
tive minimum is found. Between stages the program checks 
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the stopping conditions on the criteria function value, 
the arguments, and the maximum number of stages allowed to 
reach the minimum, and will stop if any of these conditions 
is met. The step size along a line of search is variable 
in the program, which allows the search steps to decrease 
near the minimum or within tight curves and to increase 
again once these curves are past. 

A more detailed description of the optimization method, 
the algorithm and several test cases can be found in Ref. 8. 
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V. DISCUSSION OF RESULTS 



A. VERIFICATION EFFORT 

In an attempt to verify that the spin equations of 
motion describe the aircraft motion and also to check out 
the entire computer program, several runs were made for an 
equilibrium flight condition of straight and level flight. 
The choice of this flight mode was made so that the results 
could be verified by hand while searching over a reduced 
number of independent variables in satisfying the same cri- 
teria function as discussed earlier. 

As mentioned earlier, thrust terms were not included in 
the general equations of motion, as previous studies indi- 
cated that its use is normally avoided, since flameouts are 
likely and serious engine damage can result. However, in 
order to simulate straight and level flight thrust was in- 
cluded with the very restrictive assumption that it was 
orientated in the negative X direction and passed through 
the center of gravity. 

Straight and level flight can be simulated in program 

o 

SPIN by specifying a value of ip equal to 90 , y snd Z both 
set to zero. By assuming that a wings— level, straight line— 
of— flight orientation can be obtained with neutral rudder 
and aileron control, the value of 6 and 6 are set equal 
to zero as well as the bank angle <{) . The value of the 
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radius vector (R) is arbitrary. With this accomplished for 
each run a value of the pitch angle (9) (in this special 
case 6 is the same as the angle of attack) is specified. 
Three independent variables (thrust required, elevator 
angle and airspeed) and their initial step sizes are speci- 
fied in the main program and the routine yields the value 
of these three variables that would be required for straight 
and level equilibrium flight. 

The results of these runs were verified by using an 
electronic calculator; the equations of motion were satis- 
fied. The numerical results of thrust, velocity and eleva- 
tor deflection angle were reasonable for that particular 
flight configuration. 

B. PROGRAM SPIN RESULTS 

Several SPIN runs were conducted for configuration B in 
order to compare the results with those given by Adams 
(Ref. 9) . The mass, dimensional characteristics and aero- 
dynamic coefficients for this aircraft configuration are 
listed in Appendix C. 

The results of the predicted spin characteristics for 
configuration B are shown in Table I where they are com- 
pared with the results obtained by Adams. The predicted 
spin characteristics developed by program SPIN, using the 
optimization algorithm EXTREM, are in good agreement with 
those predicted by Adams. The slight difference in the 
first, second or third decimal digit in some cases is prob- 
ably a result of the particular computer being used. 
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Adams used a CDC machine v/ith a bigger bit word than the 
IBM— 360 and therefore more significant figures could be 
handled in the numerous calculations performed. As a re- 
sult of this the minimum value of the criteria function 

—24 

obtained with program SPIN was 10 . Adams has been able 

to drive his criteria function residuals (also a sum of 
the squares but of a different set of terms) to a value of 
10-^^ 

The third grouping of results in Table I show very 
little similarity, despite the same control settings being 
used. The first steady spin predicted is an entirely dif- 
ferent equilibrium condition as can be seen by comparing 
the values of the radius (R) and angular rate (y) . The 
second spin predicted would at first appear to be the same 
as that given by Adams. The difference in values of the 
angle of attack (a) , inertial velocity (V) and radius (R) 
make it an entirely different spin condition and not just 
a poor representation of the same spin. This statement is 
justified by the closer correlation of results, in the 

other groupings, for the same value of the criteria func— 

—24 

tion (10 ) . 

The computer prediction runs for configuration A are 
shown in Table II. There were no other results to compare 
with this configuration. The spin prediction scheme was 
applied to .this particular aircraft to assure that the 
method would at least work for another aircraft. The choice 
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for this second configuration, over other possibilities, 
was made easy because the aerodynamic derivatives were al- 
ready in tabulated form as given in Ref. 7. 

The effects of density variation on the equilibrium 
spin conditions is shown in Table III. Grantham and Graf- 
ton (Ref. 8) concluded from their studies that a decrease 
in density gave a faster rotation rate (y) / higher veloc- 
ities and little change in angle of attack and sideslip. 
The parameters predicted by program SPIN verify this be- 
havior with the added effect that the radius of the spin 
tends to decrease, i.e., the spin becomes tighter with a 
decrease in density. The change in the orientation angle 
0 would indicate that the spin would be somewhat flatter 
with a decrease in density. 
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TABLE III. Effects of Density Variation 
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VI. CONCLUSIONS AND RECOMMENDATIONS 



As a result of developing and presenting the analytic 
spin prediction technique, the following conclusions and 
observations are made. 

1. The original developed equations of motion, as 
presented in Appendix A, are too highly coupled 

to be an effective tool in this type of prediction 
scheme. In this light, those equations can be con- 
sidered, at most, as an academic exercise leading 
to a rather unique formulation. 

2. The comparison of results demonstrated the accu- 
racy of the method and the effects of density 
variation study showed one possible use of the pro- 
gram. The capability of this method to predict 
steady spin conditions is limited only by the accu- 
racy with which the mass, physical characteristics 
and aerodynamic coefficients are determined. 

3. Due to the direct modeling in a cylindrical coordi- 
nate reference frame this method has a certain 
advantage over that demonstrated by Adams (Ref. 7). 
The number of equality constraints on the solution 
required by this prediction method is far fewer, 
which allows for a much simpler computer formula- 
tion and equally accurate results. However, this 
method uses more computer time than Adams ' in 
converging to the equilibrium spin condition. 

4. As with any optimization technique the importance 
of equally weighing the residuals of the equations 
is of utmost importance. The value of the residuals 
should be as close as possible to each other, in 
order of magnitude, for proper convergence of the 
optimization algorithm. Unequal weighting of any 

of these residuals cause these terms to dominate 
and influence the search in a direction that doesn't 
actually contain the desired solution minimum. Many 
hours of research time were lost by not recognizing 
this important fact. 

In light of the demonstrated usefulness of the analytic 
prediction scheme the following follow-on work could prove 
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valuable in better understanding the spin problem. It is 
recommended that: 

1. The use of rotary— balance aerodynamic data be in- 
vestigated as a possible way of more closely relat- 
ing the aerodynamic coefficients of the tunnel 
model to that experienced by the real airplane. 

The use of this type of data will allow at— desk 
calculation of the steady spin conditions. This 
particular spin estimation method, if used to pre- 
dict the initial conditions for program SPIN, 
would reduce any sensitivity to the guessed initial 
conditions and also reduce the required computer 
convergence time. This particular method of spin 
estimation is given in Ref. 9. 

2. The equations of motion could be linearized in 
order to obtain some information about the stability 
of the predicted equilibrium spin condition. 

3. This prediction scheme could be coupled with exist- 
ing fixed— base spin simulators to provide the 
starting initial conditions for the spin simulation. 
With this the recovery characteristics could be 
investigated, recovery techniques perfected and the 
necessary pilot interactions evaluated and studied. 
This one tool alone could provide valuable insight 
into the problem areas and shows much promise in 
future work. 
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APPENDIX A 



FULL EQUATIONS OF MOTION 



The following six equations completely describe the 
rigid body motions of an aircraft, using body axes, as 
formulated in a cylindrical coordinate reference frame. 



Z EQUATION 



"^^o = “l3^X “23^Y "^33^Z 



R EQUATION 



mR - m^^R = + a22^Y “32^Z^ 

Y EQUATION. 

mR(Ry + 2Ry) + Ix^P“l3 ^‘^13^ ^Y ^*^^23 ^'^23^ 

• • • 

+ I {— p(6cos(t)sin0 + 4)sin(f>cos9 ) + 
aZ 

• • • _ 

p cos(f)cos6 — 9rcos0 — rsinO} = 

^^“ll^X “21^Y “31^Z^ “l3^^X “23^Y ^33^Z 

Ip EQUATION 

^X^^“l3 ^“l3^ ^^“23^ 

l 2 (ra^^ + ra^^) + I^^ ^“P ( + 

(J>sin({)cos0 ) + p cos(f)cos0 — 0r cos0 — r sin0} = 



“l 3 ^X *^ 23 ^^Y *^ 33 ^Z 
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0 EQUATION 



I (q coS(}) — q(J)Sin(j)) — I (rsincj) + r<J)cos<{)) 

+ (ti»+Y) {l„pcos9 + I qsin(J)sinO + I rcosij)sin6} 

A Y Li 

- 1^2 sin())p) + ('i'+Y) 

+ rcos0 (ij;+Y) ) = cos<}) — M sincf* * 

Y Li 

(p EQUATION 

Ij^P + ~ p(«i)Cos<}) + sin(j)cos0 (ii»+Y) ) ) 

= M 

X 

where , 

= cosipcosd 

= — ij)siru|)cos0 — 9cosi/;sin0 
a ^2 ~ sim|;cos0 
a ^2 “ i/)cos4;cos0 — 6sini/)sin0 
“ — sin0 
“ “0COS0 

Q 22 _ “ cosij;sin0sin(j) — sini/)cos<}) 

• • • 

tt 22 ^ = — i/»sini|;sin0sin()) — 0cosij;cos0 sincj) 

• • • 

— (})cosij;sin0cosi|; — i/;cosi|;cos(}) + c})siru|;sin(J) 

tt 22 “ sini|;sin0sin(})+ cosi/;cos<j) 

• • • 

^22 ~ 0sini|)cos0sin(j) 

• • • 

— (})sini|;sin0cos(}) — i|;siruJ;coS(|) — (})cosi|)Sin(|) 

^23 ~ 

• • • 

*^23 ~ <j)cos0cosi|) — 0sin0sin(}) 

= cosij;sin0cos(j) + sini|)sin(j) 

• • • 

“ — i|;sini|;sin0cos(}) + 0cosi|;cos0cos(|) 

• • * 

— (j)cosij'sin0 sintf) + t|icos!{)sin(i) + (+'5imI;cos<}i 
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a ^2 “ sini|^sin6cos(}) — cosij;sini(j) 

a ^2 ~ '['COSijjsinOcosii; + 0sintj;cos0cos(}) — c{)sin(})sin0sin(!) 
+ i|;sint/jsin(J) — cjjcosipcoscJ) 

= cos0cos(j) 

• • • 

a ^2 ~ ~ 0sin0cos(J) — cj)sin(()cos0 

p = ^ — sin0 (t|)+Y) 

p = (}) — ij;0cos0 — i|;sin0 — y0cos0 — ysin0 

• . • * 

q = 0cos(}) + cos0sin({) (ip+y) 

• ** •• •• 

q = — 0(})sin(j) + 0cos(j) — i{)0sin0sin(}) — i|;4)cos0cos4) 

+ t|;cos0sin(}) — y0sincj)sin0 + y^cos0coscj) 

+ ycos0sin({) 

r = COS0 cos(}) (t|;+y ) — 0sincf) 

• • • • • 

r = — i|)0sin0cos(j) — i|;(j)cos0sin(}) + tj;cos0cos(}) 

• • • • • • 

— 0(j)COS(}) — 0sincj) — y0sin0cos({) 

— y^cos0sinc() + ycos0cos4> 
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APPENDIX B 



\ 

PROGRAM SPIN 



A. INTRODUCTION 

The intent of this appendix is to 'serve as a users 
manual or guide to the use of Program SPIN. An alphabeti- 
cal listing of the computer variables and their meanings 
is provided. A brief description of the main purpose of 
each subprogram as well as an extensive explanation on pre- 
paring the data deck is included. A listing of the program 
with numerous comments is included along with several func- 
tional flow charts to help in understanding the structure 
and relationships of the subroutines. A sample output is 
given and the storage and time requirements of the program 
are discussed.. 



B. LIST OF COMPUTER VARIABLES 

The following- is an alphabetical listing of the compu- 
ter variables and their meaning. 

All— A33 Direction cosines 

A( ) Working matrix for Subroutine EXTREM 

ALPHAD Aircraft angle of attack, degrees 

ALPHAR Aircraft angle of attack, radians 

B Wing span, ft 

BETAD Angle of sideslip, degrees 

BETAR Angle of sideslip, radians 
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GEAR 


Wing chord, ft 


CL 


Total roll moment coefficient 


CM 


Total pitch moment coefficient 


CN 


Total yaw moment coefficient 


COEFA( ) 


Aerodynamic coefficients which are a function 
of angle of attack only 


COEFAB( ) 


Aerodynamic coefficients which are a function 
of angle of attack and angle of sideslip 


CPHI 


Cosine of PHI 


CPSI 


Cosine of PSI 


CR( ) 


Working vector of interpolated coefficients 


CTIIETA 


Cosine of THETA 


CX 


Total X force coefficient 


CY 


Total Y force coefficient 


CZ 


Total Z force coefficient 


D 


1. O-(ZXI) **2/XI*ZI 


DDELA 


Incremental step size of aileron control, de- 
grees 


DDELE 


Incremental step size of elevator control, 
degrees 


DDELR 


Incremental step size of rudder control, degrees 


DELTAA 


Aileron control deflection, degrees (Positive 
with right aileron trailing edge down) 


DELTAE 


Elevator control deflection, degrees (Positive 
with trailing edge down) 


DELTAR 


Rudder control deflection, degrees (Positive 
with trailing edge left when viewed from above) 


DFMAX 


Stopping condition of the functional variation 
in Subroutine EXTREM 


DGAMD 


Incremental step size of GAMDOT, rad/sec 
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DPMI 


Incremental step size of PHI, degrees 


DPSI 


Incremental step size of PSI, degrees 


DR 


Incremental step size of R, ft 


DRDOT 


Incremental step size of ROOT, ft/sec 


DTHETA 


Incremental step size of THETA, degrees 


DTHRUS 


Incremental step size of THRUST, lbs 


DXMAX 


Control parameter for Subroutine EXTREM 


DX( ) 


Array of step sizes of the independent variables 


DZDOT 


Incremental step size of ZDOT, ft/sec 


F 


Function to be minimized 


FMAX 


Desired minimum value of the optimization 
function 


FOPT 


Optimum value of the function being minimized 


G 


Gravitational acceleration, ft/sec2 


GAMDOT 


Rate of change of the cylindrical orientation 
variable y, rads/sec 


HRHOS 


0. 5*RH0*S 


IW 


Control parameter of Subroutine EXTREM for 

writing instructions 

+1 All output suppressed 

+2 Final output only 

+3 Output at the end of each stage 


K 


The number of independent variables being 
searched over 


LMAX 


Stopping condition on the number of stages of 
Subroutine EXTREM (A negative sign indicates 
that a minimum is sought) 


NPRINT 


Print indicator, set in MAIN and passed to 
Subroutine DATAIN to echo check the aerodynamic 
coefficients ( =0, No print; =1, Prints) 


P 


Angular velocity about the X body axis, rad/sec 


PEAR 


P*B/2V 
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PHI 


Euler orientation angle (p , degrees 


PHIR 


Euler orientation angle p, radians 


PSI 


Euler orientation angle p, degrees 


PSIR 


Euler orientation angle p, radians 


Q 


Angular velocity about the Y body axis, rad/sec 


QBAR 


Q*CBAR/2*V 


R 


Cylindrical radius, ft 


REAR 


R*B/2*V 


RDOT 


Rate of change of the cylindrical radius, 
ft/sec . 


RES( ) 


Residual of any of the spin equations 


RHO 


2 4 

Density of air, lbs— sec /ft 


RZ 


Angular velocity about the Z body axis, rad/sec 


S 


2 

Wing area, ft 


SPHI 


Sine of PHI 


SPSI 


Sine of PSI 


SSAM 


Aircraft mass, slugs 


STHETA 


Sine of THETA 


SUMRES 


Sum of the squares of the individual equation 
residuals 


THETA 


Euler orientation angle 6 , degrees 


THETAR 


Euler orientation angle 0, radians 


V 


Velocity of the vehicle center of mass, ft/sec 


VSQ 


Velocity squared 


VXB 


Vehicle velocity component along the X body 
axis, ft/sec 


VYB 


Vehicle velocity component along the Y body 
axis, ft/sec 
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VYVB 


VYB/V 


VZB 


Vehicle velocity component along the Z body 
axis, ft/sec 


VZXB 


VZB/VXB 


X 


Working vector that contains the independent 
variables 


XFORCE 


Aerodynamic force along the X body axis, lbs 


XI 


Moment of inertia about the Z body axis, 
slugs— ft^ 


XMOM 


Aerodynamic moment about the X body axis, 
ft— lbs 


YFORCE 


Aerodynamic force along the Y body axis, lbs 


YI 


Moment of inertia about the Y body axis, 
slugs— ft2 


YMOM 


Aerodynamic moment about the Y body axis, 
ft-lbs 


ZDOT 


Vehicle velocity in the cylindrical axis 
direction, ft/sec 


ZFORCE 


Aerodynamic force along the Z body axis, lbs 


zi 


Moment of inertia about the Z body axis, 
slugs— ft^ 


ZMOM 


Aerodynamic moment about the Z body axis, 
ft-lbs 


ZXI 


Cross product of inertia about the X and Z 
body axes, slugs— ft^ 



C. PROGRAM DESCRIPTION 

Program SPIN consists of a main controlling program and 
six subprograms as illustrated schematically in Figure Bl. 

The MAIN program's most important function is to set 
the independent search variables and their initial step size 
Otherv/ise the MAIN program is basically an input/output 
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routine. The initial conditions and optimization control 
parameters are read into common memory. The initial con- 
ditions are printed out as are various final conditions 
of the airplane. 

Subroutine DATAIN is called by the MAIN program to read 
in the aerodynamic coefficients. Coefficients that are a 
function of both alpha and beta are read into the COEFAB 
array and those that are a function of alpha only are read 
into the COEFA array. DATAIN will write out all of the 
coefficients if the variable NPRINT is set equal to one. 

Subroutine EXTREM is the optimization program which 
determines the minimum of the criteria function without 
having to calculate any derivatives. The single most im- 
portant function of this subroutine is in determining the 
optimum direction of search since the criteria function 
being minimized resides, and is evaluated, in Subroutine 
FN. The returned values of the criteria function are used 
to interpolate a parabola along a line of search, to 
determine the extremum on this parabola and to determine 
new directions of search. 

Subroutine FN evaluates the criteria function value 
through calls to Subroutine CONST, Subroutine CLOKUP and 
Subroutine SSEQNS. In addition, boundaries on the search 
are -imposed due to the limited range of angle of attack 
and sideslip angle aerodynamic coefficients available. 

Subroutine CONST calculates various equation parameter 
values needed in the equations of motion. Values of 
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angle of attack and sideslip are determined so that the 
limits of the available coefficients aren't exceeded. 

Subroutine CLOKUP takes the passed values of angle of 
attack and sideslip and performs a linear interpolation 
on all the aerodynamic coefficients. 

Subroutine SSEQNS contains the six’ equations of motion 
of the airplane equated to the residual (RES(I), 1=1,6) 
vector. The residuals are evaluated and passed, through 
common memory, to Subroutine FN so that the criteria func- 
tion can be evaluated. 

Functional flow charts of the MAIN program and each 
subroutine are shown in Figure B2 through B8 . 
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Figure Bl. Schematic of Program SPIN 
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Figure B2. MAIN Program Functional Flow Chart. 
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Figure B3. Subroutine DATAIN Function Flow Chart. 
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Figure B4 . Subroutine EXTREM Functional Flov; Chart. 
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Figure B4 . (Continued) 
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Figure B5. Subroutine FN Functional Flow Chart. 
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START 







CALCULATE SPIN 
EQN. VARIABLES 







RETURN 



Figure B6 . Subroutine CLOKUP Functional Flow Chart 




Figure B7. Subroutine CLOKUP Functional Flow Chart 
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Figure B8 . Subroutine SSEQNS Functional Flow Chart. 

D. PREPARATION OF THE DATA DECK 

The following is a guide for preparation of the data 
deck used by computer program SPIN. 

CARDS 1 and 2 — The initial conditions on GAMDOT (rad/ 
sec), R(ft) , RDOT (ft/sec) , PHI (deg), THETA (deg), PSI (deg) , 
and DELTAA(deg) on card one and DELTAR(deg), DELTAE(deg), 
THRUST (lbs) and ZDOT (ft/sec) on the second. Data is input 
as floating point according to the format ( 7F11 . 5/4F11 . 5 ) . 

CARD 3 — The controlling parameters for Subroutine 
EXTREM. DFMAX, DXMAX, FMAX, LMAX and IW are read in ac- 
cording to the format ( 3E10 . 2 , 15 , 12 ) . 

CARD 4 — Aircraft identification name (e.g., CONFIGURA- 
TION A) starting in column one. 

CARDS 5 and 6 — The mass and physical characteristics 
of the aircraft, atmospheric density and gravitational 
acceleration values. Card 5 contains: SSAM (mass, slugs). 
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2 — 

S (wing area, ft ) , B (wing span, ft), CBAR (c, ft) and G 

2 

(gravitational acceleration, ft/sec ). Card 6 contains: 

2 2 2 
XI (I , slugs— ft ) , YI (I , slugs— ft ) , ZI (I , slugs— ft ) , 
X y z 

2 

ZXI (I , slugs— ft ) and RHO (atmospheric density, slugs/ 

X z 

3 

ft ) . Data on each card are input according to format 
(5F16. 8) . 

CARDS 7 — 412 — (alternating three card series) . These 
cards contain the tabulated aerodynamic coefficients that 
are dependent on both angle of attack and sideslip. The 
first card in the series is as follows: in columns 1&2 

the coefficient number (I) , columns 3&4 blank, column 5 
the index on beta (IB) , columns 6—10 blank and the remain- 
ing 70 columns, in 10 column increments, are the coeffi- 
cients indexed on alpha. The second card of the series is 
divided into eight 10 column increments and the third card 
is divided into four 10 column increments of aerodynamic 
coefficients indexed on alpha. Alpha is indexed 19 times 
for each beta index which corresponds to a range of alpha 
from 0 to +90 degrees, in 5 degree intervals. Beta is 
indexed nine times which corresponds to a range of beta 
from —40 to +40 degrees, in 10 degree intervals. 

The first set of three cards, within a particular coef- 
ficient grouping, corresponds to an angle of sideslip equal 
to —40 degrees, where the first coefficient in columns 
11—20 corresponds to an angle of attack equal to 0 degrees, 
columns 21—30 corresponds to an angle of attack equal to 
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5 degrees and so on in 5 degree increments of alpha until 
90 degrees is reached in columns 41—50 of card three. 

The second set of three cards contains coefficients over 
the alpha range for beta equal to —30 degrees. The remain- 
ing series have coefficients for beta equal to —20°, —10°, 
0°, +10°, +20°, +30° and +40°. The next grouping of 27 
cards (nine three card series) contains another aero- 
dynamic coefficient over the range of alpha and beta and 
so on in 27 card groupings until all 15 alpha and beta 
dependent coefficients have been prepared. Format for the 
three card series is ( 12 , 2X, II , 5X , 7F10 . 7/8F10 . 7/4F10 . 7 ) . 

The input sequencing of groupings is very important as 
the machine computation is dependent upon the assignment 
of a particular tier in the three-dimensional COEFAB array 
to the proper aerodynamic coefficient. Table Bl gives the 
proper assignment. 

CARDS 413—429 — (alternating three card series) . These 
cards contain the tabulated aerodynamic coefficients that 
are dependent on alpha only. The first card in the series 
is as follows: columns 1&2 the coefficient number (I), 

columns 3—10 blank and the remaining 70 columns, in 10 col- 
umn increments, are the coefficients indexed on alpha. The 
second and third cards of the series are the same as de- 
scribed above for cards 6—411. Again alpha is indexed 19 
times corresponding to a range of alpha from 0 to 90° in 
5° increments. Columns 11—20 of card one corresponds to 
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the coefficient for alpha equal to zero degrees, columns 
21—30 for alpha equal to 5°, and so on in 5° increments 
until alpha equal to 90° is reached in columns 41—50 of 
card three. Format for these three card series is (12, 8X, 
7F10 . 7/8F10 . 7/4F10 . 7) . As before, the input sequence of 
the coefficient groupings is very important for proper 
interpretation by the program. The assignments are given 
in Table Bl. 
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TABLE Bl. Coefficient Assignment 



Coefficient 


storage 

Location 


Coefficient 
Print Out 


Interpolated 

Location 


^1 


COEFAB 


( 


1 


CR( 1) 


^m 


COEFAB 


( 2,_,_) 


2 


CR( 2) 


C 

n 


COEFAB 


( 3,_,_) 


3 


CR( 3) 


C 

y 


COEFAB 


( 4,_7_) 


4 


CR( 4) 


c 

X 


COEFAB 


( 5,_,_) 


5 


CR( 5) 


C 

z 


COEFAB 


( 6,_,J 


6 


CR( 6) 


6r 


COEFAB 


( 7,_,J 


7 


CR( 7) 




COEFAB 


( 8,_,_) 


8 


CR( 8) 


^n6r 


COEFAB 


( 9,_,J 


9 


CR( 9) 


^y 6a 


COEFAB 


(10,_,_) 


10 


CR(IO) 


*^16a 


COEFAB 


(!!,_,_) 


11 


CR(ll) 


^n6a 


COEFAB 


(12,_, J 


12 


CR(12) 


c . 

z 6e 


COEFAB 


(13,_,J 


13 


CR(13) 


C 

m6e 


COEFAB 


1 

'l 

1 — 1 


14 


CR(14) 


c . 

x6e 


COEFAB 




15 


CR(15) 


C 


COEFA 


( 1,_) 


16 


CR(16) 


P 












COEFA 


( 2,_) 


17 


CR(17) 


P 










c 

n 


COEFA 


( 3,_) 


18 


CR(18) 


P 










C 


COEFA 


( 4,_) 


19 


CR(19) 




COEFA 


( 5,_) 


20 


CR(20) 


r 










C 

n 


COEFA 


( 6,_) 


21 


CR(21) 


r 










C 

z 


COEFA 


( 7,_) 


22 


CR(22) 


q 










c 

m 


COEFA 


( 8,_) 


23 


CR(23) 


q 










c 

X 


COEFA 


( 9,_) 


24 


CR(24) 


q 
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ooo or>ooo ooooo ooooo oooo oooo oooooooooooo 



E. PR0GR7VM LISTING 

jJ: A :;c j}: V V ^ Xc 2^; 3}s 2{c sjc 3 ^ r^c ^ :}c :^c a}c a|x aj: ajc a^ jJ: a^: aj: :& ^>: :(; a';; a|: 

^'/SIK PROGRAM =(= 

=(: 

PLRPCSE SERVES TO SET THE INCEPENCENT VARIAELES TC 

BE SEARCHED OVER AND CCKTRCLS THE READING CF =<'■ 
THE INITIAL CONOITICNS AND CPTIHIZATICN =!= 

CONTROL PARAMETERS. ALL STRUCTURE FOR THE 
PRINTED OUTPUT, '.v I T H THE EXCEPTION CF THE * 

OPTIMIZATION RESULTS, CCRES FROM THIS 
PROGRAM. =!= 



IMPLICIT REALMS! A-H,0-Z ) 

EXTERNAL FN 

DIMENS ION X( 11) ,DX( ID ,A( 11, 1A> 

CCKMCN/ SPIN/ GAM DOT .GAMMA, ROOT , P H I , THET A , P S I .CELTAA ,0EL 
ITAE, DEL TAR, THRUST , V.COEFAB (15,20,10),COEFA(20,9),SSAM, 
2S ,6,CBAR,XI , YI , ZI , ZXI , G ,RHO ,A I l,AU,AL3,A2l,A22,A23,A3 
31 ,A32 ,A35 ,VSQ, VXl , VY 1, VZl ,P , C , R Z , P EAR , QB A R , RBAR.HRHCS, 
4CR (24) , RES (6) , AL PHAO , A L PHAR , bcT AD , EET AR , S F E T A , C BE T A , S A 
5 LPHA,CALPHA,SGAMMA,CGAMMA, SPSI , CPSI , STh E T A , C THET A , SPHI 
6 , CPHI , R .XFCRCE ,XFORl , Y FCRC E , Y FC R 1 , ZFORC E , Z F C R 1 , XMQM, YM 
7 0 M , ZMO M , Z 0 QT , NP R I N T , N AM E ( 4 ) 



DATA INITIAL STEP SIZES OF THE INDEPENDENT 
VARIABLES 



DATA DGAM/C. l/.DPSI/ .Oi/,DZDQT/5 .0/ ,DPHI/ .C 1/, OR/5 .C/, 
IDTHETA/ .31/,DDELE/C.l/,DTHRUS/LCC.C/,DRDCT/5.0/ 



READ THE INITIAL CCNDITICNS ON THE INCEFENDENT 
VARIABLES AND OPTIMIZATION CONTROL PARAMETERS 



PEAD(5 , 100) GAMOCT, R.RDOT , PM , T H ETA , PS I , 0 E LT A A , DELTAR , 
IDELTAE, THRLST ,ZDOT 
REAC(5, 101) DFHAX, CXMAX,FMAX,LMAX, IH 

SET PRINT INDICATOR FOR THE AERODYNAMIC 
COEFFICIENTS: NPRINT=0, SUPPRESS PRINTING 

NPRINT=l, PRINT CCEFFICIENTS 



NPR1NT=0 

SET THE NUMBER QP INDEPENDENT VARIABLES (K=) AND 
ASSIGN EACH NAME TO THE X ARRAY WITH THE CORRES- 
PCNOING STEP SIZE TO THE DX ARRAY. 



K = 6 

X( 1)=GAMD0T • 

X(2)=R 

X(5)=ZDCT 

X(4)=PSI 

X(5)=THETA 

X(6)=PHI 

DX( 1) =DGAM 

CX(2) =DR 

DX(3) =OZOCT 

DX(4)=0PSI 

CX(5)=DTHETA 

DX(6)=DPHI 

CALL TC SUBROUTINE CATAIN TC BEAD MASS, PHYSICAL 
CHARACTERISTICS AND AERCCYNAMIC COEFFICIENTS 
OF THE AIRCRAFT. 

CALL DATAIN 

WRITE OUT THE INITIAL CCNDITICNS 

WRITE! 6, DC) NAME 
rtFITE (6 ,102) 
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c 

c 

c 

c 

c 

c 

c 



=«■ MAIN PROGRAM CCNTINLED * 

WRITE(6,103) 

WPI7£(6,10A) GAMDCT ,R,ROGT ,PHI , THETA 
WFIT5(o,l05) 

URIT£(6,10A) PSI ,DELTAA,OELTAR,CELTAE,THRLST 
WRITEI6, 106) 

WPITE(6,119) ZOOT, PHC , OHM AX , CXM A X , FMAX 
UR ITE ( 6 , 107) 

WPITE(6, 10c) LMAX, lU 
UPITEI 6 ,109) 

CALL THE OPTIMIZATION ALGORITHM TO FINE THE 
STEADY SPIN PARAMETERS. 

CALL EXTREMI FN,K, X, CX , A , OFMAX , D XM A X , FMA X , L M A X , FOP T , I U) 

WRITE CUT THE FINAL STEADY SPIN PARAMETERS. 

WRITE(6,110) NAME 
WRITE(6,111) 

W P I T E ( 6 , 1 1 2 ) 

WRITE (6 ,10A ) ALPHAC, 6ETAD, GAMDOT , THETA, PH I 
WRITE! 6, 113) 

WRITE(6 ,10A) PSI , GAMMA,R, THRUST , RHC 
WRIT£(6,114) 

WPIT£(6,1C4) V,VX1,VY1,VZ1,P 
WRIT£(o ,115) 

WRITE(6,116) Q,RZ,RDGT 

WRITE(6,117) XFORCE, YFCRCE, ZFORC E , XMQM , YM C M ,ZMCM 
DC 2 1=1,0 
WRITEI6, 120) 
wRirt(6, lie) 

STEP 

FORMAT (7F 11. 5/4F1 1 .5) 

FCfMAT ‘ ‘ 

FORMAT 
FCRM “ 

1 • FHI , . 

2MFT/SEC) .w. 

104 FORMAT (T54,5F11.5) 

1C5 FORMAT ( / / , T5 8 , ' P S I ' , T69 , • D ELT A A ' , T30 , ■ L» t L i A 

1 • CELT Af: • , T 10 2, • THRLST ',/T57,'(DEG)' ,T69 , ' (DE 
2* (DEG) • ,79 1, • (DEG) • ,T 102, ' (LBS ) • , / ) 

106 FORMAT ( // , T 5 8 , ' ZD 0 T • , T69 , ' RH'O • , T 80 , ' DF M A X ' , 

I'CXMAX' ,T1C2,*FMAX',/T55,'(FT/SEC)',/) 

107 FORMAT (// ,T5S, 'LMAX' ,T69, ' IW* , / ) 
ice FORMAT (T57, I 5,T68, 13,/) 

109 FORMAT (5(/),T67,‘ - EXTREM ITERATIONS ',//) 

lie FORMAT ( ' 1' ,8(/) ,T73 ,4A4,// ) 

111 FORMAT (T67-‘- - - pI'Vai dao 

112 format (//,. 

1 ' theta* ,T1 02 

2 * (PAD/SEC) ■ 

112 FORMAT (// 



I ,RES ( I ) 
FOPT 



ICC 

101 

1C2 

102 



lAT (7F 11. 5/4F1 1 .5) 

'AT (3E10.2, 15, 12) 

'AT (T67,*- - - INITIAL CONDITIONS - - -') 
1AT (// ,T58, 'GAMDOT* ,T71, ‘R* , T eC , * ROOT * , T9 1 , 
■ * ,T 102, ‘THcTA * ,/T56 , * (RAD/SbC ) * ,To9 , ' (LEG) * 
r/SEC) * ,T90, * (DEG) * ,T102 ,* (DEG) * , / ) 

'AT (T54,5F11.5) 



,T78 



* DELTAR * ,T9 1 , 
■(DEG)*,T80, 

rT91, 



IT 1 C2 , *KnC 



rLKnAI l///fJDOt 5L/r (Jr Iht oWiJPI 

L*OF THE RESIDUALS' ,£17.5) 

FORMAT (T53,2F11 .5,2X,J£11 .4) 
FORMAT (//,T6 9, *RES( * , I 1, * ) = * , E li 
F^r. 
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Q * :C: Xc ^ ^ t ^ ^ ^ a; A al' 4c ^ 4^ >1- 5 ;: a? 3^1 it ^ >{= « 5{: 3^ 



SLBRCUTINt CATAIN =1-' 

PLRPCSE READS THE PASS AND CIMENSICNAL CHARACTERISTICS 
OF THE AIRCRAFT AS KELL AS THE AERCCYNAMIC ^ 
COEFFICIENTS. 






.J. ^ «|L« •X* O. <U. OU .t. ’J* w- 

*lr •'.' ^ 'T* -T' '»• -r 



«l. «f . U. a. U, xU a. «l. iJ. %L. . .1 .J. U, .J, O. ^ .Jf ^ 1.)^ ^ Of Ou «L. 

nr -r -r ^ -n* T* V V n* ^ ^ ^ "r T* 



SLBRCUTINE DATAIN 
IPFLICIT REAL=<^8( A-F,3-Z) 

CCFNCN/SPIN/GAMDCT ,GAPPA ,RDGT , P F I , TF ET A , P S I .CELT A A , CEL 
IT A E .DEL TAR , IHRUST , V ,C QE F AB ( 15 , 2 C , 1 C ) , C OE F A ( 2C , 9 ) , S S A P , 
2S , L-.CbAR ,X I ,Y I , ZI , ZX I , G , RHO , A L 1 ,• A 1 4 . A 1 3 , A 2 1 , A 2 2 . A 2 5 , A3 
2I,A32,A33,VS;',VXl,VYL,VZl,P.C,RZ,FtAR,»^EAR,RBAk,FRFCS, 
4CR ( 2A ) , RES (6i , AL PH AC , A L PHAR , c ET A L , 0 ET AR , S E E TA , CdE T A , S A 
5LPFA,CALPhA,SGAPMA, CGAPPA, SPS I .CPSI.STHEIA.CIFETA.SPHI 
6.CPH1 , R ,XFCF C£ , XF CR 1 , Y FCRCE , YFGR i , ZFORGE , ZFCRl .XF'CP, YN 
7LM, ZMCF , ZDCT , KPRINT ,NAME( A J 

READ IK THE GC KF I GURAT I ON NAPE, MASS AKC PHYSICAL 
PARAMETERS OF TFE AIRPLANE, AIR DENSITY AND 
GRAVITATIONAL ACCLERA.TICK VALUE. 



RE AD (5, 100) NAME 
REAC(5,1C1) SSAM,S,E,C8AR,G 
READ(5,101) XI ,YI ,ZI ,ZXI ,RH0 

READ TFE AERODYNAMIC COEFFICIENTS THAT ARE A 
FUNCTION OF ANGLE-OF-ATT ACK AND SIDESLIP 

CC 1 J=l,135 

1 READ(5,1C2) I , IB , ( CCEFAB ( I , lA, I B) , IA=l , 19 ) 

READ TFE AEROCYNAMIC COEFFICIENTS THAT ARE A 
FUNCTION OF ANGL E -OF-A T T AC K ONLY 

CC 2 J=l,6 

2 READ(5,103J I , { COE F A ( I A , I) , I A = 1 , 19 ) 

CC 3 J=l,3 

3 REAC(5,1C4) I , ( COEFA( I A, I ) , I A = 1 , 19 ) 

CC 4 1=1,15 

CC 4 J=1 , 10 
CCEFAB ( I ,2C, J) =0.0 

4 CCEFAB! I ,J , 10)=0 .0 
CC 5 K=l,9 

£ CCEFA( 2C,K)=0.0 

CHECK PRINT INDICATOR FOR ECFG CHECK CF DATA 

IF{NPR INT .EC.O) GO TO 9 
CC 7 1=1,15 
l\RITE(6 , 105) NAME 
K F I T E ( 6 , 1 1 C ) 

KRIT£(fc,lC6) I 
KRITE! c, 107) 

KRITE (6 , 108) 

KRITE(6,1C9) 

N = C 

CC 6 IA=1, 19 

KFITE(6,111) N,( CCEFAB (I,IA,IB) ,IE = 1,9) 

6 N = N + 5 

7 CCNTINUE 
N = C 

KRITEI6.105) NAME 
KRITEIe ,113) 

KRITcIc ,112) 

CC 8 1=1,19 

KFITE(6,111) N',(CCEFA(I,J),J=1,9) 

£ N=N+5 
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a[s :{ij5< * yr :J: jjc Jt ^ -r ?? ?:t 5}! sit ri; :?c :*: :^c 2;s 3*-. a{: * jJc :f: ilc :{c :J: : 

* SUBkOUTINe CATAIN CONTIKUEC ^-' 

S RETURN 

ICC FCPRAT (4AA) 

FCRMAT (5F16.8) 

FCF.iMAT (I2,2X,I1,5X,7F 1C . 7 / OF 1 0 . 7 /AF 1 0 . 7 ) 

FCP.PAT ( 12 ,8X, 7F10 .7/8F10 .7/AFlO .7 ) 

FCPMAT ( 12, IX,7F 11 .7/7FL1.7/5F 1 1.7) 

FCRN.AT ( • 1 • ,9(/ ) ,T63 ,4AA) 

FCRMAT ( T65, ' COEFFICI ENT *,I2,//) 

FCRMAT 



101 

102 

1C3 

ICA 

1C5 

106 

107 



1 '- 10 ' 

2 '+ 40 ' ) 
108 FCRMAT 
FCRMAT 
FCRMAT 
FCRMAT 
FCRMAT 
I'CGEF 



(T29, ' 



,T77 



O' 



BETA 

,T35, 



,T4C, 
' + 10 ' 



40' 



1C9 
1 1C 
1 1 1 
112 



(T27, ' (DEGREES ) ' 
(T29, • ALPHA' ) 

(/) 

(T3C, 12 ,4X ,9F9 .5 
(T29,' ALPHA' ,T38 
18 • , T65 , 'COEF 19* 



I — 

,T94, 

) 

, / ) 

,'CGEF 
,T74 , 



,T49, 

'+2C 



'-30' ,T58 , '- 
,T103, • 43C ' 



20' , T67 
T112, 



16 ' 
'COEF 



,T47,'CCEF 17',T56 



2C ,T32 , ‘ 



V. 1_ t: r 



2T92,'C0EF 22 • , T 1 0 1 , ' COEF 23 ' , T 1 1C , 'COEF 24',/T27 
3 ' ( DEGREES ) ' ,/ ) 

113 FCRMAT (/) 

END 



21 
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^ ^ :<< 2[^ >> : 



^ « X*' : 



X^ X'- 2^ a^ X' X: >{• 5\ X^ X: ajt :Jt Xt 3^ >",1 X< X- 



SIBRCUTINE EXTREM 



!{: 



PLRPCSE 






THE OPTIf'IZ/MICR ALGCRITRR CALL EC FFCR THE *■ 
P'AIN PROGRAM THAT RETURNS THE MINI MM FCUi\0 ’S^ 
FOP A PARTICULAR STEADY SPIN SEARCH . THE <• 
SUBROUTINE IS PROVIDED LITE AN INITIAL LINE 
OF SEARCH BY THE USER ThPCUbH THE CALLING 
ARGUMENT LIST AND THE ALGCRITHM CET ERMINES *• 
THE MINIMUM ALONG THIS LINE bY I N T E R PC L AT I NG 
A PARABOLA THROUGH THiREE CRITEPxIA FUNCTION 
VALUE POINTS DETERMINED EY CALLS TC SLcROUTINE 
FN. secondary DIRECTIONS OF SEARCH ARE 
DETERMINED BY A GRAM-SCHMIET ORT HCGCNAL IZAT ION 
PROCESS AND CONTINUES, '^ITh VARYING STEP 
SIZES, UNTIL A MINIMUM IS FOUND. * 

t- 

Wb. «A, JU S.V. «Ju 'O. O. ^ *A« ^ •.(> vl. «v ^ <Jb V O. '*# 'V' '.v «v V> 

'Z' •'r -r n' -n ^ 'T' ^ "r ^ *r* 'c 't* -*r nr ■.-r rr nr "r* -'n nr n' nr n* V ^ ^ nr ' rr n' ^ -r 



SLBRCLITIKE EX TREM ( F , K , X , DX , S , DF M AX , CXMAX , F M AX , UMAX , 
IFCPT, IK ) 

IMPLICIT REAL=i=8(A-H,0-Z) 

DIMENSION X( 1 ) ,DX( 1 ) ,S III , 1 } 

ASSIGN COUNTER AND FLAG VALUES 



L = C 
LI = 1 
N = C 

DC 1 1=1, K 
S ( 1 , 1 ) = X ( I ) 

1 S(i, 2 ) = x(n-Dxm 

CALL 1C SUBROUTINE FN TO CHECK IF THE INITIAL 
CONDITIONS EXCEED THE BOUNDARIES OF ANGLE-QF- 
ATTACK AND SIDESLIP 

CALL F(X,F2,LI,N) 

FE = F2 

IF(LI .GT.l )WRITE(6,2) 

2 FCRMATI IX ,' INITIAL ARGUMENTS OUTSIDE BDUNCARIES') 

3 I F (KC .GE.K .OR .KC .LT.O.OR.L . EQ. 0 ) KC = 0 
KC=KC+1 

S( 1,3I=C. 

CC A I = I, K 

S( I ,4) =S( I ,1)-S ( I ,2 ) 

A S( i,3) = S( 1 ,3)+S( I ,A)--^=f=2 
S ( l,3) = CSGRT(S( 1 , 3 ) ) 

IrUABS(lW).GE.3)KP ITE(b,23)L,N,S(l,3),F2, 

1 ( 1,X(I ) , I ,CX( 1} , 1 = l,K) 

CHECK THE STOPPING CONDITIONS 

IF(L.GE.IA6S( LMaX ) .OR .5 (1 , 3 ) .LT . CXMAX .OR . DA £S ( FF-FGPT ) 
1 .LT .DFMAX.AND. L. GT .O.CR.ILI.GT.l.AND.L.EC.O.CR.FOPT. 
2LE .FMiAX .AND. L .GT .0 )GO TC 22 
IF (K.EC .U GCTC9 

COMPUTE THE ORTHOGONAL SECONDARY DIRECTIONS BY A 
GRAM-SChMICT C RT HOGON AL I ZAT I CN PROCESS 



CC 8 J=2,K 

KC=-2 + J-*KC 
IF (KD.GT.K )KD=KD-K 
S( J,3)=C. 

DC 7 I =1 , K 
S( I ,J+5)=0. 

IF ( I .EC .KC)S( I , J + 3) = S( 1,3) 
JM= J-1 

CC 6 JK=1,JW 
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SLifcROUTINE EXTREM COMIRUEC 

>■>■■* «»•- 'Jf ^ >'* -^o ^ »JU V' *Ax •!*» >f «■’' '5!r •»'' 'V ''■* ^ **» Oi Ju «*, ^ %).. ^ OU W>t 

^ •T “T •T "r ^ ^ n' 'T* " ' *r* V - -r* ^ *r* -r -r 't* '(^ t' 'O ^ *T* '<• •*r^ -r *x ~ /* "ir* €> ^ ^ ^ ',•*»» »r ' t- > 'T n* ^ 



6 S ( I , J + 3 )--=S( I , J+3 )-S (KD, JK+3 l,2)/S(JK,3)=!‘S(I,JK + 3) 
1/S( JK,3) 

7 S(J,3) = S( J,3)+S( I , J + 3)-=*2 
S( J,5J=CSQPT(S( J,3 ) ] 
iF(S( J,3) .LT.1.0-3 0)GOTC3 

e CCMINUE 
9 DC 10 1=1, K 

IC S( I ,2]=S( I ,1J 
L = L + 1 
FF=FGPT 



STAGE CO LOOP 

DC 21 M=1,K 
DC I I I = I , K 

1 1 S (I ,M+3)=S( I,y+3)/S(M,3)*0X(K) 

IF( IW.GT.O )LI=3 

12 IF ( IW .GT .0 )L I = LI-l 
IJ=LI 

DC 13 1 = 1, K 

X(I) = S( I, U-S(I,,M + 2) 

13 S(I ,K+3)=S( I,l)-X( n 

CALL SUBROUTINE FN TO DETERMINE THE CRITERIA 
FUNCTION VALUE AT X-DX 

CALL FIX, FI, LI, N) 

EC= U 

lA DC 15 1=1, K 

X(I)=S( I,1HS(I ,N + 3)/BC 
15 SI I ,M+3I =XI i )-SI I , 1) 

IF IDA8S leO) .GT . 1 . I ) GOTO 20 

CALL SLBRCUTINE FN TO DETERMINE THE CRITERIA 
FUNCTION VALUE AT X+DX 

CALL FIX,F3,LJ,N) 

IF CUTSIDE THE BOUNDARIES GO BACK AND ADJUST 
THE STEP SIZE 

IF ILI + LJ.EG.^» )G0TC12 
IF ILJ .GT .2 )BC=-A . 

IFILI .GT .2 )E*0 =+4 . 

IFILI.GT .2.0R.LJ.GT.2)GCT014 
IE ST=C. 

INTERPOLATE A PARABOLA TFROUGF THE THREE 
FUNCTIONAL VALUE POINTS: X-CX, X, X+CX 

DC 18 1=1, K 
XI I )=SI 1,1) 

IFIDABSISI I,M+3) ) . LT . I .D-30 ) GOTO IE 
S I I ,M+3 ) = S I I ,M+3 )/LI 

1FIDABSI2.^^F2-F1-F3).LT.1.0-30)GCTC18 
XIl)=StI,l) + SII,M + 3)/DABSIF 1- 2 . F2 + F3 ) ~ I F 3- F 1 ) / 
11SIGNI2 ,LNAX) 
le ST = ST + i XI I )-SI 1,1) )"-*2 

1 F I 16 .-ST .LT .DXIM )**2 )DXIM) =DXI N )/A 

I F I ST .LT .4CC.*DXIM )*^-2 .AND .DABS I 2. *F2-F 1-F3 ) .GE. I . C-30 
.1)GC TO 20 
CC IQ 1”1 K 

IFIDABSISI I,M+3) ) . LT . 1 . C-30 ) GOTO 19 



INTERPOLATE FOR THE APPROXIMATE MINIMUM POINT 

X I 1 ) = S I I , 1 )+OS IGN I S I I , iM+3) , IF3-F 1 ) /SI I , M + 3 ) ) ^ 

1 IS ICN (2C , UMAX ) 

19 CONTINUE 



73 



:*i :{« ;ft 3): 3^: ij 2*: ^ tjc 5V sj'. it 3l: :V :^^ 5“r 2V * 5> y- 3^ 3fe r^? 34 A 5.V aj: :<« H ii i: i:i' 

^ SUEROUTINE EXTREM CCNTIKUEC 

3{l:^r5}:3>3^^|c2''^3:. 5iV>,V5;-*3:t5j«3^. 4t:*;y:3;i3,V3i?3^F3;c3;s3V5!c3;c;{?:^c:^:;::it3;s5*<3;:5X:iC3!:^:<:3^3!<3{i*3;cy^3{:3^:^c:{*y^:J:2;x:^>;c3;<^ 

CX (N) = DX(M )*2 . 

2C LI^+1 
EC=-BO 

iF( CABS (EG ) .GT . 1 . 1 )DX (M) = OX(K )/3 . 

C 

C DETERMINE THE MINIMUM VALUE 

C 

CALL F(X,FCPT,LI ,N) 

IF(LI.G7 .1 )LI=10 
C 

C SET THE NEXT CPTHOGONAL SEAPCF DIRECT ICN AND 

C STEP SIZE 

C 

I F( iSIGN< 1 ,LMAX)>'<( FCPT-F2I .LT.-CAES(FE-F2)^A. .AND .L I . 
INE.10JLI=2 

IF(LI .GT.l .AND.DABS(BO) .GT.1.DGCTG14 
I F ( LI .GT . i IGOTOlb 
FE = F2 
F2=FOPT 
DC 21 1=1, K 
2 1 S ( I , L ) = X( I ) 

GCTC3 

22 IF(IABS(Ils).EC.2)WPITE(6,23)L,N,S(l,3),F2, 
l( I ,X( I ) , I , CX( I ) , 1= l,K ) 

23 FCFMAT (// ,T57, ' STAGE N'C . ' , I 5 , 20X , ‘ TR lAL N C . • , I 6 , / / T 57 , 
1 'CL=' , E 15 . 5, //T57, 'FUNCTION VALUE = • , E 15 . 5 , / // ( T 57 , 
2'AB(‘,I2,‘)=',E15.5,5X,'DS(‘,I2,')=',E15.5)) 

RETURN 

END 
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^ ^ 5^ X* ^ ^ a ^ 2 !* 5';: si: ^ j!; * 5^ 4c 5}: 5}C aj: ^ a*: 3^ =»c 5^ 

SLdHCLTINE FN * 

•ju 

FLRFCSE O^LLEC BY SLEFCUTINE EXTREF' FOR C E T E RR I NAT I CN^ 
CF THE CRITERIA FUhXTICN VALUE. SLERCUTIKE - 
FN MAKES CALLS TO SLbRCUTIKE CONST, CLGKUP, T= 
AND SSEQKS TO CETERMNE TFE CRITEhU FUNCTION^ 
VALUE. A CHECK TO INSURE THAT THE VALUES CF * 
ANGLE-OF-ATT /5CK AND SIDESLIP IS PRCVICED SC 
THAT A SEARCH tiEYDNC THE CCEFFICIENT TABLE v 
LIMITS IS PROHieiTEC. IN SUCH A CASE, CCNTPCL 
IS RETURNED TO SUBROUTINE EXTREM SC TEAT TFE * 
CORRECT STEP SIZE CAN BE ADJUSTED TC VALUES * 
LYING WITHIN TFE TABLE. ^ 

JL. 

J: Xi ^ 4 3^ V ^ >1^ * ^ X: s}; 4 : s? 4 4 4 V X' ❖ 4 4 X* 4 4 X* 4- 4 X' 4 4 4 4 X: 



SLBRGUTINE FN(X,F,LI,N) 

IMPLICIT REAL«6( A-H,C-Z) 

CIMENSICN X( L) 

CCMMDN/SPI N/GAMDCT ,GAMMA, RDOT , PF I , THETA, PS I ,DELTAA ,CEL 
IT AE,DE L TA R , THRUST , V ,CD£FAa (15,20,1C),COEFA(2C,9),SSAM, 
2S , B,CBAR,XI , YI , ZI , ZXI ,G,RHO ,A I l,AU,A13,A2i,A22,A22,A3 
31 ,A32,A3^,VSQ,VX1 , VY1,VZ1 ,P,C,R2 ,PEAR,Q8AR,RBAR,HRFCS, 
ACRl 2Y) ,RES (6) , AL PH AD , ALRHAR ,BETAD,BETAR,SEETA,CBETA,SA 
5LFFA,CALPFA,SGAXMA,CGAMFA,SPSI , CPSI,STHETA,CThETA, SPHI 
6,CPFI,R,XFCRCE, XFC R i , Y FCRC E , Y FOR 1 , ZFCRCE , 2 F C R1 , XMC F , YM 
7GM,ZMCM,ZDCT,NPRINT ,NAME(4) 

GAMDOT = X( 1 ) 

R=X( 2) 

^CCT=X(3) 

FSI=X(4) 

TFETA=X(5) 

PF I=X(6 ) 

CHECK TO SEE IF THE BOUNDARIES ARE VICLATEC, IF 
SO, CONTROL IS RETURNED TO SUBROUTINE E>TREM WHERE 
THE STEP SIZE IS ADJUSTED CF NEW ARGUFENTS 
DETERMINED THAT DC NOT VIOLATE THE BCLNLARIES 

CALL CCNST 

I F( ALPHAO.lt . 0.0 .OR. ALPHAD.GT.90 . .CR.BET AD . LT .-40 . .CR. 
lEETAD.GT .40. JL I = L I + 1 
IFILI .GT.l JRETURN 

CLOKUP MAKES TFE INTERPCLATED COEFFICIENT VALUES 
AVAILABLE IN CCMMCN MiEMCPY 

CALL CLGKUP 

SSECNS MAKES THE RESIDUALS AVAILABLE IF CCMFICN 
MEMORY 

CALL SSEQNS 

CRITERIA FUNCTICN D ET E R F I NAT I C N 

F = RES( 1 J^RESI 1)+RES(2)*RES( 2 ) +R E S ( 2 ) >!^RE S ( 3 ) +R E S ( 4 J * 
IRES (4J + RES (SJ^IRESISJ+RESISJ^J-RES (fc ) 

N = N+1 

RETURN 

END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



3^ 

SL6PCLTINE CCNST * 

:U 

PURPCSE DETERP'INES VARICUS ECL'ATIGh PARAMETERS NEECEC* 
IN EVALUATING THE EOUATICNS OF NCTICN. VALLES 
OF ANGLE-CF-ATT ACK AN'C SIEESLIP ARE PASSEE * 
TO SUEROUTINE FN FOR A LIMIT CHECK ANC ALSC 
TO SUEROUTINE CLOKUP FOR CCEFFICIENT INTER- * 
RELATION. 



5\ ^ V 5 ; 5*; ::c * :ic 5^: J'i sf: 5J. 3tv ift :J5 >1' 5,‘: ::ic 5k 5,r 5^ 5 ;-. v : 



: 5|: 5 ;- 5^ 5^ 5^: :^s s}-, 5)? 5|c 5{, :J; 



SLEROUTINE CONST 
II'PLICIT REAL*E( A-F ,0-Z ) 

CCMMCN/SFIN/GAMUCT ,GAPMA ,RDOT , PF I , TFET A , F S I .DELTA A, CEL 
IT AE, DEL TAR .THRUST , V .COEFAB ( 15,2C,1CJ,COEFA(20,9J,SSAR, 
2S,e,CBAR,XI.YI,ZI,ZXI,G,RHO,Ali,A12,A13,A21,A22,A23,A3 
21,A32,A23,VSJ,VXl,VYl,VZl,P,G,R2,PtAR,UdAR,PBAR,FRFOS, 
ACR(24) ,RES(t>) , AL P h AD , A LPHAR , BET A C , EET Ak , S E E T A , CBE T A , S A 
5LFFA.C ALPFA.SGAMMA.CGAMPA, SPS I.CPSI.STHETA.CTFETA.SPFI 
6 ,CPHI , R .XFCRCE .XFCRl.YFCRCE.YFORl.ZFORCE.ZFCRl ,XMCP , YM 
7CM ,Zf-'CF , ZDCT, NPR INT ,NAME(4 ) 

ThETAR=Th ETA/57 .295779 
PF IR=PHl/57. 295779 
PS IR=PS 1/57.295779 
STFETA=CSIN(THETAR ) 

CTFETA=DCCS( THETAR ) 

SPFI=DS JN(PHIR) 

CFFI = DCCS( PHIF.) 

SFS I = DS IN( PSIR) 

CFSI=DCCS( FSIRJ 
A13=-STFETA 
A23=SPhp'^CTHETA 
A32=CPHi-CTHETA 
A12 = CTFETA''^SPSI 

A22=SPh I>^-STFETA^'SPS I + CPh I-CPSI 
A32 = CPHl-STh£TA^SPSI-SPFI='^CPSI 
A 1 UCTHETA^i-CPSI 

A21=SPHI-STHETA-CPSI-CPHI*SPSI 
A 3 1 =C PH I * S THE TA*C PSH-SPHI^-SPSI 
HRFOS= .5^^RHQ'^S 

VSG = (GAMDQT»RJ=^-(GAMC0T*R) +ZCCT*ZCCT + RDOT=^RCCT 
V=OSQRT ( VSG) 

P = A 13=^-GAM0GT 
e=A23«GAMCCT 
RZ=A33-GAMDCT 
PEAR=0 . 5=^P^B/V 
CeAR = 0. 5vC^=CB AR/V 
REAR = C. 5-RZ='-d/V 
ARGUM=GAMCCT^P/V 
GAMMAR=CARCGS ( ARGUM ) 

GAMMA=GAMMAR*57.295779 

VXe = Al lvGAMDOT=i^R + A 12 R CGT + A 13*Z CGT 

VYB=A21*GAMDOT=!=R + A22'-i=RCGT + A23«ZCGT 

VZD = A31='-GAMOOT*P+A3 2'>=RDGT + A23vZCGT 

VX1=AI l-VXB+A21-VYe + A31'-VZ6 

VY1=AL2*VXB+A22Y VYe+A32-VZB 

VZ l = AI3^VXE+A23^VYS + A33=^VZB 

VZXe=VZB/VXB 

VYVB= VYB/V 

ALPhAR=DAT AN( V2X3) 

EETAR = OARS INIVYVB ) 

ALPFiAD = ALPHAR*57.2 957 79 
BETAD=BETAR"57. 295779 
SALPHA=CS!N( ALPriAR ) 

CALPEA=DCGS(ALPhAR) 

SEETA = DSIMcETAR) 

CEETA = CCGS (BETAR J 

RETURN 

ENC 
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StePCLTINE CLCKUP * 

* 

PLRPCSE TAKES THE PASSED VALLES CF ANGLE-C F-ATTACK =!= 
AND SIDESLIP AND PERFORMS A TWD-F,AV LINEAR * 
INTERFOLAT ICN CF ALL IS CCEFA6 (ALFFA AND * 
BETA DEPENDENT COEFFICIENTS ) ARRAY CERIV- 
ATIVES. A LINEAR INTERPOLATION IS PERFORMED =:• 
CN THE 9 CCEFA ( ALPHA CNLY DEPENDENT COEF- 
FICIEMTS ) ARRAY DERIVATIVES. THE INTERPG- 
DATED COEFFICIENT VALUES ( CR(I),I=i,LA ) ARE* 
STORED IN CCMNCN MEMORY FOP USE BY SUBROUTINE* 
SSEONS IN EVALUATING THE ECUATICNS CF MGTICN.* 






SIBROUTINE CLCKUP 
IMPLICIT REAL='--8( A-h,0-Z) 

CCKMGN/SPIN/GAMDGT ,GAMMA, RDCTtPPI , THETA ,PSI,DELTAA,DEL 
1TAE,D£LTAR,THRUST, VtCOEFAB ( 15,20,101, CO EFA(2C, 9) , SSAM, 
2S,B,C6AR,Xi,Yi,ZI,Z>:i,G,PhC,Ail,A12,A13,A21,A22,A2-2,A3 
3 1,A32, A23tVSC,VXl,VY 1 , VZ 1 , P , G , R Z , P BAR , Q B AR , P B A R , H RhC S , 
AC F (24) , RES (6) , ALPFAC , ALPHAR , B ET A D , B ET AR , S E F T A , CBET A , S A 
SLPHArC ALPHA,SGAHHA ,CGAMMA,SFSI , C FS I , STHE 1 A , CTHET A , S PF I 
6, CPHi , K,XF0HCE,XFCR1, YFGRCE ,YFOR 1, ZFORCE , ZFCR1,XMGM, ,YM 
7CM,ZMCM ,ZDCT , NP R I N T , N AM E ( 4 ) 

AI=G.2*ALPFAD+1.G 
EI = 0.i’^BETAD + 5.0 
IA1=AI 
IB 1 = BI 
I A2 = IAl4l 
IE2=IBl+l 
DX=bI-IEL 
CY=AI-I Al 

INTERPOLATE ANGLE -OF-ATT ACK ANC SIDESLIP 
DEPENDENT COEFFICIENTS 



DC 1 K= I, le 

C1=C0EFAB(K, lAl, IB1) + CX*(C0EFAB(K,IA1, IB2)- 
ICCEFABI K , I AL , IBI ) ) 

C2 = C0EFAB (K, I A2, IBI) +0X* { C 0 E F AO ( K , IA2, I B2 )- 
1CCEFAB( K , I A2, IBI ) ) 

1 CR (K)=C 1 + 0Y-(C2-C1 I 

INTERPOLATE ANGLE-OF-ATT ACK DEPENDENT CECFFICIENTS 

CC 2 M.= l,9 ■ 

I =M+15 

2 CR( I)=CCEFA( IAI,M)+DY* (COEF A (I A2 , M ) -C OE F A ( I A 1 , M ) ) 
RETURN 

END 
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SLBRCUTINE SSECKS * 

PLPPCSE CALCULATES THE TOTAL FORCE AND MCNEAT COEF- « 

FICIENTS ANC THE RESIDUALS ( RES(I),I=1,6 ) * 

CF THE SIX EOLATIONS CF ACTION. THE RESIDUALS 
ARE MADE AVAILABLE IN CCAACN MEMCRT FCR - 

SCBRCUTIME FN TO EVALUATE THE CRITERIA =!- 

FUNCTION. * 






SUBROUTINE SSEQNS 
lAPLICIT REAL*8(A-H,0-Z) 

CCAA.CN/SPIN/GAMDCTrGAMAA,RDOT ,PH I , THETA , P S I iDE LT A A , CEL 
ITAEtDELTAR ,THRUST, V,COEFAB (15,20,1C),COEFA(2C,9),SSAA, 
2S ,B ,CDAR ,X I ,Y I » 21 i 2X1 , G ,RhO , A1 1 ,A12,A13,A21,A2 2,A22,A3 
2l,A32,A33,VSQ,VXl,VVI,VZI,P,C,R2,FEAR,0BAR,FdAR,HRHCS, 
4CP (2A ) , RES (6 ) , ALPHAD, ALPHAR , BET AC , EETAR , SBETA,CBETA , SA 
5LPHA,CALPHA,SGAAPA,CGAAAA,SFSI , C F S I , ST H ET A , CTHETA, SPHI 
6 tCPHI , R,XFORCE, XFCR L, YFCRCE , YFOR I , 2 FORCE , 2FCR 1 , XACA , YM 
7CA,2MCA,2DCT, NPRINT,NAME(4) 

CALCULATE THE TOTAL FORCE ANC AOPENT CCEFFICIENTS 

CX = CR(4 ) + CR( 15)=!--DELTAE + CR{24)=5=QEAR 
CY^CR(4)+CR{ L6J''’^PBAR + CR( 19 ) »:^R BAR-J-CR ( 10)*CELTAA + 

ICR (7)=s=OELTAk 

C2 = CR(o H CR( 13)^=DELTAE + CR(22)=!=CeAR 
CL=CR( L )+CR( 1 1 ) ^D£LTAA + CR (8 ) =^DELTAR+CR( 17 )*PBAR + 
1CR(2C)>:^RBAR 

CP = CR (2 )+CRI L4)*DELTAE+CR( 23)*QBAP 
CN = CR (3 ) + CR( 12)'-^DELTAA+CR(9 )^DELTAR+CR( 18 )>»FBAR + 
lCP(21)vRbAR 
C=1.0-ZX i’^ZXI/IXI-'i'ZI ) 



CALCULATE THE AERCCYNAMIC FCRCES AND FCRENTS 

XFCBCE:=HRHCS=^VSCYCX+THRUST 
YFCRCE^HRHGS*VSw*CY 
ZFCRCE=- HRHOS^fVSO-CZ 
XPCP==HRHCS'=VSG-B*CL 
Y M C M = H R H 0 S V S C C b A R C M 
2 P C P = H P H C S * V S Q - B * C N 

XFCRl = All-XFC)KCEi A2l^YFGRCE + A3I=f 2FCRCE 
YFCR1 = A 12*XFORCE + A22*YFGRCE + A32=>ZFCRCE 
2FCRl = A13’"'XFORCE + A23=^YFORCE + A33’='ZFCRCE+SSAP*G 

CALCULATE THE EQUATION RESIDUALS 

R£S( l)=XFORI/<SSAP*V) 

RES(2 } = ( YFOR 1-SS AM'i^G A PCCT*6 AMDOT R ) / ( SS A P =» V ) 
RtS(3)---2HaRl/ (SSAP^V) 

RES(4) = (({YI-Zn/XI-(2XI*ZXI)/(XI=«=ZI))-:'Q*PZ4(l.-(YI-XI 
1 )/2I 2X I^=P”Q/XI H ( HRHOS*VSQ-B/Xl )=t(CL+(ZXI*CN/ZI))) 
2/C 

RES( 5 ) = HRHOS=i-'VSQ^CBAR»!^CM/YI + ( ZI-XI )'-'P*RZ/Y I + (ZXI/Y I )=«= 

I (R2'-=RZ-F=^P ) 

RES<c.) = nZXI*ZXI/(XI"‘ZI)-(YI-XI)/ZI)*P-G+((YI-ZI)/XI- 

II . )-ZXi*Q4K^/ZI + ( HRHQS*VSg=^e/ZI ) '"M ZX P^^C L /X I -fCN ) ) / D 
RETURN 

END 
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SAMPLE OUTPUT 



F, 



CONFIGURATION A 



INITIAL CCNOITICNS 



GANCOT 

(RAC/SEC) 

1.01000 

PSI 
( CEG ) 

75.CC000 



R 

(FT ) 

23.CC000 

CELTAA 

(DEG) 

15 .OOCOO 



ROOT 
(FT/ScC ) 

0.0 

DELTAR 

(OEG) 

-30.00000 



PHI 

(DEG) 

5.00000 

OELTAE 

(OEG) 

-25 .00000 



TFETA 
( DEG) 

-35 .OOCOO 

TFRUST 
( LBS) 

0 .0 



ZDCT 

(FT/SEC) 

340.00000 



RHG OFPAX DXMAX FMAX 

C.C0089 0. 10000-60 C.lCOOC-60 0.1000D-6C 



LRAX Ik 

-1050 2 



- EXTRER ITERATICNS - - - 



STAGE NO. 656 TRIAL NG. 11830 



DL= 0.0 



FUNCTION 


\/ALU£= 0 


.25 


726C-23 






AP ( 


1 ) = 


C .S6C60D 


OC 


0S( 


1) = 


0.2220AD-16 


AP ( 


2) = 


C .69878D 


C2 


DS( 


2 ) = 


0.5A2 iCD-2l 


AR ( 


3) = 


C .34345D 


C3 


DS( 


3 ) = 


C.338BLC-19 


AR ( 


4) = 


0 .20299D 


02 


DS( 


A ) = 


C. I lD-15 


AR ( 


5) = 


-0.39 L63D 


02 


OS ( 


5 ) = 


0 . 1 L3O90- 14 


AR ( 


6) = 


0.31776D 


C2 


DS( 


6 ) = 


0. S47C3D-20 
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CONF IGCRAT ICN A 





FINAL 


PARAMETEFS 


— — — 




ALPPAO 

(DEG) 


BET AC 
(CEG) 


GAMCCT 

(RAD/SEC) 


THETA 

(CEG) 


Phi 

(DEG) 


4£. 38319 


1C .45719 


0.98060 


-39 .16334 


20.29875 


PS I 
(CEG) 


GAMMA 

(CEG) 


R 

(FT) 


THRUST 
( LBS ) 


RHO 


69 .87808 


84 .81609 


31.77595 


0.0 


C.CCC89 


V 

( FT/SEC} 


VXl 

(FT/SEC) 


VV 1 

(FT/SEC > 


VZl 

(FT/SEC) 


P 

(RAD/SEC) 


34A .86437 


31.15947 


-c. coccc 


343 .45381 


0.61928 


e 

( RAC/SEC) 


RZ 

( PAC/SEC) 


RCCT 
(FT/SEC ) 






C .26376 


C .71309 


o.c 









XFORCE= 
VFORCE= 
2FCRCE= 
XMOM= 
YFCM = 
ZMCM= 


2992.653C0 
-7896.07153 
-68433.28272 
5437 .88A27 
-127702.55131 
42513.27919 




RES( 1 ) = 


-0.760C5D- 12 




PES(2)= 


-0.12141C-12 




RES(3) = 


0. 13899C-1 1 




FES(4)= 


-0. 322260-13 




FES (5) = 


0 .38337C-13 




RES (6 > = 


0.I1154C-12 


SLR CF 
OF THE 


THE SCLARES 
RES ICLALS 


0 .25726C-23 
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G. PROGRAMMING CONSIDERATIONS 



Program SPIN requires 66K bytes of storage. When buffer 
space is taken into account, a request for 80K bytes should 
be made when submitting the job for processing. 

The time requirement is variable and is very dependent 
on the closeness of the initial conditions to the steady 
spin condition. Convergence has been obtained with a total 
CPU time of a little over one minute to as high as 15 min- 
utes; this time includes about 25 seconds for the compile 
and link steps. 

The optimization scheme loop time is very rapid. The 
time required to complete one stage is 0.2 second, which 
requires 18 criteria function evaluations. The criteria 
function evaluation occurs at the rate of one per .0111 
second and this accomplishes, as Figure Bl indicates, the 
evaluation of the equation constants, interpolation of the 
24 aerodynamic coefficients, and the evaluation of the 
equations of motion. By looking at the listing for each 
of the subroutines involved, one can better appreciate the 
speed with which the program operates. 

The variable time requirement of the program can create 
a problem to the researcher. More often than not, the 
turn around time for any program is determined by the time 
and storage requested. As a general guideline to the use 
of program SPIN, the following is suggested: 

a. Set LMAX (the maximum number of stages) to a value 
of 1,050. 
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b. Request 80K storage and a four minute run time. 

c. Check results of the run. If the value of the 

—5 

criteria function is greater than 10 , guess 

another set of initial conditions and rerun. Re- 
sults indicate that if the value of the criteria 
function hasn't dropped below the above value in 
the four minutes alloted that an increase of time 
will not yield a steady spin condition. The pro- 
gram in this case is converging on a local mini- 
mum which is not indicative of a steady spin. If 

the value of the criteria function is less than 

—5 

10 , rerun the program with LMAX set at 2,500 

and a requested time of nine minutes if greater 
accuracy is desired. 
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APPENDIX C 



AIRCRAFT DATA 

The following aircraft data was used in the study of 
steady spin conditions. Table Cl tabulates the mass and 
dimensional characteristics of the two aircraft and the 
following pages list the aerodynamic coefficients for each 
configuration . 

TABLE Cl. Mass and Dimensional Characteristics 



Characteristics 


Configuration A 


Configuration B 


m (slugs) 


1554.00 


771.81 


s (ft^) 


525.00 


695.00 


b (ft) 


63.00 


38.00 


c (ft) 


9.04 


23.76 


I (slugs-ft^) 


53100.00 


13600.00 


I^ (slugs-ft^) 


299000.00 


128000.00 


I^ (slugs— ft^ ) 


338750.00 


138000.00 




12480.00 


4340.00 
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